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Abstract 


A numerical method based on b-spline polynomials was developed to study incom- 
pressible flows in cylindrical geometries. A b-spline method has the advantages of 
possessing spectral accuracy and the flexibility of standard finite element methods. 
Using this method it was possible to ensure regularity of the solution near the origin, 
i.e. smoothness and boundedness. Because b-splines have compact support, it is also 
possible to remove b-splines near the center to alleviate the constraint placed on the 
time step by an overly fine grid. Using the natural periodicity in the azimuthal direc- 
tion and approximating the streamwise direction as periodic, so-called time evolving 
flow, greatly reduced the cost and complexity of the computations. 

A direct numerical simulation of pipe flow was carried out using the method 
described above at a Reynolds number of 5600 based on diameter and bulk velocity. 
General knowledge of pipe flow and the availability of experimental measurements 
make pipe flow the ideal test case with which to validate the numerical method. 
Results indicated that high flatness levels of the radial component of velocity in the 
near wall region are physical; regions of high radial velocity were detected and appear 
to be related to high speed streaks in the boundary layer. Budgets of Reynolds stress 
transport equations showed close similarity with those of channel flow. However 
contrary to channel flow, the log layer of pipe flow is not homogeneous for the present 
Reynolds number. A topological method based on a classification of the invariants 
of the velocity gradient tensor was used. Plotting iso-surfaces of the discriminant of 
the invariants proved to be a good method for identifying vortical eddies in the flow 
field. 
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Chapter 1 


Introduction 


This chapter presents an overview of numerical investigations of flows in cylindrical 
geometries, especially pips flow. Some background is also included on the important 
features of pipe flow. The objectives of the present study are laid out, and we conclude 
with an outline of the present work. 


1.1 Numerical method 

1.1.1 Background 

Since the early seventies, numerical simulations of fluid flows have become an in- 
valuable tool in the study of turbulence. Direct numerical simulations (DNS), or the 
simulation of fluid flows without the use of any model to resolve the smallest scales 
of motion, are the tool of choice when the detailed physics of fully turbulent flows are 
of primary importance. 

The expertise and the computational resources have grown so much in the last 
few years that numerical simulations do more than just supplement experiments, 
but are now used as an independent means of exploring the physics of more and 
more complex flows. Flows, most of them incompressible, with one inhomogeneous 
direction are not uncommon; the homogeneous directions are assumed to be periodic 
which greatly reduces the complexity of the computations. Both unbounded, or free 
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shear, and bounded flows have been successfully computed using DNS; among more 
recent work, Moser and Moin [25] conducted a DNS of a curved channel, Moin and 
Kim [24] and Kim, Moin and Moser [17] a straight channel, and Rogers and Moser [38] 
a fully turbulent self-similar mixing layer. O’Sullivan and Breuer [30, 31] studied the 
growth of disturbances in pipe flow, and Sondergaard [42], studied the effect of initial 
conditions on the growth of a plane wake. These are typical of what are essentially 
two types of direct numerical simulations: the first is a transitional simulation where 
special care is applied to the initial conditions and where transition to turbulence is 
the ultimate goal of the study; the second, and more common, is a simulation of a 
fully developed turbulent flow where only the end-state is desired, with little or no 
attention paid to the transient state. 

Yet, because of the cost involved in carrying DNS, practical engineering applica- 
tions are still out of reach (and will be for decades). Because the smallest scales gets 
smaller as the Reynolds number increases, only lower Reynolds number flows in sim- 
ple geometries are attainable with present day computers. Still, this does not imply 
that DNS do not have any practical use; as an example of usefulness of DNS, results 
from such simulations are often used to “fine-tune” turbulence models used in other 
types of simulations commonly used in day to day engineering (see for example Rodi 
and Mansour [36] or Mansour Kim and Moin [22]), such as large eddy simulations 
(LES) (e.g. Mansour, Ferziger and Reynolds [21]) and mostly, Reynolds averaged 
Navier-Stokes (RANS) equations. A subgrid scale model is used to resolve the finer 
scales of motion in the case of LES, and a turbulence model is used to resolve most 
of the unsteadiness in the case of RANS. 

1.1.2 Survey of numerical methods 

Spectral methods are used to carry out DNS because of their high degree of accu- 
racy, giving them the ability to faithfully capture all scales of motion; derivatives 
of flow parameters, such as vorticity, are also represented very accurately. These 
methods work especially well on smooth and continuous fields, typical of incompress- 
ible turbulent flows. In the homogeneous directions periodic boundary conditions are 
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applied and Fourier transforms are used to simplify the equations. In addition, the 
use of fast Fourier transforms (FFT) which make the implementation very efficient. 
However, when walls are present, or more generally when a direction is not periodic, 
Fourier transforms are no longer applicable. One must then resort to other kinds of 
expansions in the cross-stream or normal to the wall direction (e.g. Spalart et al. 
[44]). 

In order to be successful, a method used to solve the incompressible Navier-Stokes 
equations in cylindrical geometries must address two principal problems: first, im- 
posing the continuity equation accurately, and second, accounting for the coordinate 
singularity which appears at the centerline (r = 0) when working in cylindrical coordi- 
nates. As we shall see, resolving problems associated with the coordinate singularity 
is the outstanding difficulty in the design of the numerical method. The problem 
of imposing the continuity equation arises when solving the incompressible Navier- 
Stokes equations, since the continuity equation appears as a kinematic constraint on 
the velocity field. Problems linked to imposing continuity constraint have been well 
documented, and the reader can consult Canuto et al. [7], Hughes [13] or Johnsson 
[15] for more on the subject. 

Several methods have been employed to study fully developed turbulent flows in 
cylindrical geometries. All these methdos use Fourier transforms in both the stream- 
wise and azimuthal directions, and some other form of expansion in the radial direc- 
tion. It should be noted that periodicity of the azimuthal direction is natural and not 
an approximation, contrary to the streamwise direction for which periodicity implies 
a flow with infinite extent. 

For their solution of transitional pipe flow, Leonard and Wray [20] (see also 
Leonard [19]), used divergence-free expansion functions in the context of a (non- 
Galerkin) weighted residual method. Using solenoidal, or divergence-free expansion 
functions alleviates the need to explicitly solve the continuity equation; thus, only 
two of the velocity components are independent. With the no-slip boundary condi- 
tion built into the expansions and using integration by parts, the pressure term drops 
out of the equations, and only two unknown velocity components remain. For the 
radial direction they used shifted Jacobi polynomials. Treatment of the coordinate 
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singularity was also built into the vector expansion functions using a priori knowl- 
edge of the behavior of the solution near the origin. This approach guarantees the 
solution to be regular * , i.e. values and derivatives are all bounded and smooth. In a 
slight variation, Moser, Moin and Leonard [26] used Tchebyshev polynomials instead 
of Jacobi polynomials to solve a fully turbulent curved channel flow. Numerically this 
flow is somewhat simpler to solve than pipe flow, since the presence of an inner wall 
eliminates the coordinate singularity. 

In analyzing the growth of both linear and nonlinear disturbances in pipe flow, 
O’Sullivan and Breuer [30, 31], and earlier, Orszag and Patera [29], used a collocation 
method with scalar expansion functions instead of vector expansions, i.e. each velocity 
component and the pressure were represented separately. Tchebyshev polynomials 
were used to represent the radial direction. Time marching was achieved using the 
splitting method of Orszag and Kells [28] where the pressure gradient and viscous 
terms were integrated in a three-step process. The nonlinear term was integrated 
using a second order Adams- Bashforth scheme. As was the case with the method of 
Leonard and Wray, it a priori knowledge of the behavior of the velocity and pressure 
near the origin was built into the expansions. 

Zhang et al. [55] also used a spectral method for their simulation of fully devel- 
oped pipe flow. Similarly to Leonard and Wray [20], they used Jacobi polynomials 
as the expansion basis although not in the context of a divergence-free formulation. 
Time marching was carried out using a fractional step method, with a mixed explicit- 
implicit second order scheme. Their method also offers the possibility of resolving 
inflow-outflow boundaries, or spatially evolving flows, although details on the treat- 
ment of those boundary conditions were sketchy and no results using them published. 

Eggels et al. [12] used a second order finite difference technique in the radial 
direction to simulate a fully developed turbulent pipe flow. Even though finite differ- 
ences are much simpler to implement than a spectral method, they induce numerical 
diffusion and dispersion, especially at high wave number (small scales). For a given 


Certain authors prefer the terms analytic or holomorphic to regular. 
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accuracy, many more grid points must be used compared to a spectral method. Ulti- 
mately this reduces the highest Reynolds number that can be adequately simulated 
due to memory constraints of the hardware. Eggels et al. [12] implemented a way to 
alleviate the strong restriction on the time step produced by a concentration or clus- 
tering of azimuthal modes at the centerline. They treated terms involving azimuthal 
derivatives implicitly, and all other terms explicitly with a leap-frog scheme for the 
advective terms, and a lagged forward Euler scheme on the viscous terms. The values 
at the origin are evaluated using a first order extrapolation; they claim that first order 
extrapolation does not degrade global accuracy because near the center, the velocity 
and pressure fields are smooth. It should also be noted that in the context of finite 
differences, ensuring that the velocity and pressure are bounded and regular near the 
origin is trivial since only their values, and not derivatives, have to be accounted for. 

Verzicco and Orlandi [49] also developed a second order finite difference scheme 
for cylindrical geometries. They treated the coordinate singularity by introducing a 
radial flux (ru r ) on a staggered grid, i.e. only the radial flux is evaluated at the origin. 
The equations were then solved with a fractional-step and approximate-factorization 
methods, yielding a second order method in both space and time. This method was 
used by Orlandi and Fatica [27] in simulating fully developed turbulent flow in a 
pipe with rotating wall, and by Verzicco and Orlandi [48] for a transitional round jet. 
However, contrary to Eggels et al ., no mention is made of any procedure to alleviate 
the constraint on the time step introduced by the increased azimuthal resolution near 
the origin. 


1.2 Pipe Flow 

Pipe flow is without a doubt the most studied flow, and the one with the broadest 
engineering applications. Reynolds himself [34] developed his famous scaling param- 
eter (c.f. the Reynolds number), by conducting experiments on pipe flow. As such, 
pipe flow represents the ideal test case on which to apply a new numerical method: 
the availability of ample experimental measurements should make validation simpler. 
However, not all issues regarding pipe flow are resolved, and the present results should 
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address those as well. The next few paragraphs give the salient features of pipe flow 
and how it contrasts with another well studied flow: channel flow. 

Pipe flow differs in many ways from channel flow (the flow between two flat plates). 
The first difference was observed by Patel and Head [32] who found that in the 
low Reynolds number range (Re), which for the purpose of the present study is the 
relevant range, fully developed pipe flow fails to conform to the standard logarithmic 
law. They found the Reynolds number had to reach a value of Re = 10000 (based on 
the bulk velocity and pipe diameter) before the mean velocity profile would match 
the log-law, versus a much lower Re = 3000 for channel flow (based on the distance 
between the plates). However, a more recent study of pipe flow by Durst, Jovanovic 
and Sender [11] using a precision laser Doppler anemometry technique (LDA) found 
that at a Reynolds number of 7442, fully developed turbulent pipe flow obeyed the 
standard log law. 

A second difference with channel flow lies in the transition from laminar to tur- 
bulent flow. Patel and Head [32] observed that pipe flow is purely laminar up to 
Re = 2000 after which it enters a transition regime usually triggered by disturbances 
at the entrance. The transition is maintained up to Re « 3000 after which the pipe 
becomes fully turbulent rather suddenly. For channel flow, they observed that the 
onset of transition appeared at Re = 1350 and the flow became fully turbulent for 
Re > 1800. The transition regime for pipe flow was studied in detail by Wygnanski 
and Champagne [52] and Wygnanski, Sokolov and Friedman [53] (see also Cantwell 
[5]). They found that the transition regime was actually a very complex process where 
two different intermittent flows were present: puffs and slugs. Puffs are present for 
2000 < Re < 2700; they are induced by large amplitude perturbations at the pipe 
inlet. Wygnanski et al. referred to the state of the flow in puffs as that of incomplete 
relaminarization. Slugs however are present at Re > 3200 and are produced by small 
perturbations in the boundary layer at the pipe inlet. The flow within a slug is similar 
to fully developed turbulence and has a well defined leading and trailing edges, and 
can extend over the entire length of the pipe. Leonard and Wray [20] and Leonard 
[19] were able to simulate puffs at Re = 2200, showing that methods designed to 
simulate time evolving flows (streamwise periodic) can capture flow features arising 
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from perturbed entry conditions. 

Lastly, pipe flow is said to be linearly stable to small disturbances. Several inves- 
tigators have focused their research on trying to find unstable linear modes in pipe 
flow (see for example Salwen et al. [40, 39], Orszag and Patera [29] and O’Sullivan 
and Breuer [30]). To this day, no unstable linear modes have been found. Channel 
flow however is linearly unstable to infinitesimal disturbances; it possesses a criti- 
cal Reynolds number of 11544 above which the flow can be unstable. For pipe flow 
the only process that sustains and enables growth of perturbations emanates from 
nonlinear interactions. 

1.3 Objectives 

The objectives of this study are threefold: 

• To develop a numerical method based on b-spline polynomials which would 
have an accuracy comparable to standard spectral methods and the flexibility 
of finite elements or finite difference methods. 

• To study the fundamental behavior of wall-bounded turbulence in a cylindrical 
geometry, such as pipe flow. 

• To establish the benchmark in pipe flow simulation. 

Also, in order to assess the present methodology, several comparisons will be made 
with various experiments and other direct numerical simulations. 

1.4 Outline 

The present work is essentially divided in two parts: a detailed explanation of the 
numerical method, and analysis of the results. 

Chapter 2 introduces the numerical method, including some basic facts about b- 
splines, the time marching method and the technique used to relax the resolution 
near the origin; it concludes with some test cases used to validate the method. 
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Chapter 3 presents the results from the pipe flow simulation; results include, mean 
flow properties, spectra, turbulent intensities, budgets of Reynolds stress transport 
equations and structure tensors. 

Chapter 4 presents a study of flow structures using a topological method based 
on a classification of the invariants of the velocity gradient tensor. 

Chapter 5 concludes with some final observations and recommendations for future 
work. 

Appendix A gives a short summary of b-splines and their properties. 

Appendix B exposes in great details the problem of the coordinate singularity and 
the steps taken to ensure smoothness and boundedness of the velocity field near the 
origin. 

Appendix C shows how the method presented in chapter 2 and appendix B was 
implemented. This appendix is especially useful if subsequent modifications should 
be made to the computer code. This appendix outlines the entire structure of the 
code. 

Appendix D gives the exact solution to Stokes problem used to validate the 
methodology and includes some background on linear stability. Those exact solu- 
tions are used in the validation section of chapter 2. 



Chapter 2 
Methodology 


In this chapter, the numerical method that was used to solve the incompressible 
Navier-Stokes equations is introduced, including some coding strategies. Also in- 
cluded are some of the test cases that were used to validate the method. 


2.1 Approach 

The numerical method was designed for any incompressible turbulent flows, that can 
be prescribed in cylindrical coordinates (e.g. pipe flow, round jet). In other words, 
the method should have the flexibility to resolve both solid and free shear boundaries. 
Since no attempt is made at resolving inflow/outflow boundaries (spatially evolving), 
the flow is assumed periodic in the streamwise direction, or so-called temporally 
evolving flow (equivalent to a flow with infinite streamwise extent). 

The present method is conceptually similar to the approaches of Leonard and 
Wray [20] and Moser, Moin and Leonard [26], in that it makes use of divergence- 
free vector expansion functions. The major difference between the present method 
and all the ones described before, is the use of basis splines polynomials (b-splines) to 
represent the radial direction, in place of Jacobi or Tchebyshev polynomials. B-splines 
are discrete polynomials, i.e. with local support on a given interval, whereas Jacobi 
and Tchebyshev polynomials have global support. Polynomials with local support 
lead to sparse matrices that can be efficiently stored and solved. The local nature of 
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b-splines enables effortless implementation of various boundary conditions, and the 
possibility of adjusting the precision near the origin to avoid an overly restrictive time 
step. Background on b-splines can be found in appendix A. 


2.2 Governing Equations 

The starting point for the simulations are the incompressible Navier-Stokes equations 
in cylindrical coordinates 


V u 


du 

~dt 


+ u • Vu + Vp 


0 

-J— V • Vu + /e : 

ne b 


( 2 . 1 ) 

(2.2) 


where u is the velocity, p is the pressure, and / is a uniform body force (pressure 
gradient) acting in the streamwise direction which drives the flow; it is constantly 
adjusted such that the mass flow, m = irR\pUb, is kept constant. Because the density 
p and the radius of the pipe R 2 are constant, the constant mass flux condition implies 
that the bulk velocity t/j, 




u z r dr 


(2.3) 


must also be constant. 


All quantities have been normalized using the radius of the domain R 2 as the 
lengthscale, and Ub as the velocity scale (see figure 2.1) 


_ _ r° _ 2 ° ^ _ t° 

r ~R2 Z= R~2 1 = jRjlJb) 

So, the Reynolds number in 2.2 is given by 


Re b = 


UJh 

V 



p = 


p° 
pU z 


(2.4) 


where v is the (constant) kinematic viscosity. 
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Using the following identity 

V(u-u) 

u • Vu = — u x u> H 

where the vorticity, u>, is defined as 

weVxu, (2.5) 

the nonlinear advective term in 2.2 can be rewritten to yield 

^ + VP = ux W + -J-V • Vu + fe z (2.6) 

at Re b 

with 

u • u 


P = p + 


2 


(2.7) 
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2.3 Numerical Method 

2.3.1 Weak Form 


Before the Navier-Stokes equations can be solved, they must be transformed in a way 
which is amenable to numerical procedures. Let v be the numerical approximation to 
u, which will consist of a truncated expansion in terms of divergence-free vector func- 
tions (i.e. V ■ v = 0), and let P be the numerical approximation to P . Furthermore, 
let £ be any vector function representable by an additional finite set of divergence-free 
vector functions (V • £ = 0). Spaces which encompass such functions can be written 
as 

V = {v | v € H 1 (G), V • v = 0 in G, v = 0 on r« 2 } (2.8) 

W = \t € H 1 (G), V • £ = 0 in G, £ = 0 on r* 2 } (2.9) 

where r fi2 represent sthe surface r = R 2 (see fig. 2.1) and H n is the Sobolev space of 
degree n. A function, say /, is said to belong to H n if it possesses square-integrable 
derivatives up to order n (Hughes [13]), i.e. 

1 dx < oo (2-10) 

With v = 0 on Tr 2 , the no-slip boundary condition, V is specific to pipe flow. 
In general, boundary conditions of any type (Dirichlet, Neumann or mixed) could 
be built into V; boundary conditions built in V are said to be imposed strongly. As 
long as Dirichlet boundary conditions are imposed strongly on functions in V (slip 
or no-slip), W will always retain the form 2.9, i.e. with £ = 0 on T/? 2 . Functions 

in W must be homogeneous at the boundary if the pressure is to drop out of the 

formulation (more on this later). 

To ensure that v and £ are truely divergence-free, a standard result from calculus 
states that if the domain G is simply connected, then V • v = 0 in G if and only if 
v = V x *, where is a vector stream function associated with v. So, if v is to 
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be continuous (C 0 ), then must at least be C 1 *. This means that second order, or 
quadratic b-splines are the lowest admissible order that can be used. 

By substituting v for u in 2.1 and 2.6 and using the standard weighted residual 
technique which consist in taking the scalar product of the equation with a weight 
function and integrating over the domain, with £ as the weight functions, the weak 
form of the Navier-Stokes equations is obtained U 

(^^) + («-VP) = («,(vxa,)) + 1 L(t(V-Vv)) (2.11) 

where u> = V x v and the inner product (■, •) is defined as 

(a, b) = f a. - bdV (2.12) 

J G 

Since V • v = 0 by construction, the continuity equation (2.1) drops out from the 
equation set. A further simplification is possible by eliminating the pressure from 
2.11. First, by rearranging the pressure integral 

/ £-VPdV= / S7-(Pt)dV- / PV-ZdV (2.13) 

Jg Jg Jg 

Then, applying the divergence theorem and using the fact that V • £ = 0 results in 

j £-VPdV = J^PZ-ndS (2.14) 

where T = T 0 U IT,, U Tr 2 and n is the unit normal vector. Using periodicity of all 
variables in the streamwise direction, any surface integral on r 0 U IT Z vanishes. With 
only the integral on r^ 2 remaining, n = e r , and since £ € W, 2.14 is identically zero. 
Hence, the weak form reduces to 

{(, §£) = («, (v x «)) + 2L («, ( V • Vv)) (2.15) 

*More formally, with # lying in some space, say V, one could write V = | $ G H 2 {G)}. 

tNote that the body force term (/e z ) was dropped from 2.11. In the code, the mass flow is held 
constant through a computational procedure described on page 130; the /e z term is never actually 
computed. 
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By enforcing 2.15 for all £ 6 VV making up a basis for the weight functions, a 
coupled set of ordinary differential equations for the coefficients in the expansion for 
v € V are obtained, which can be solved using standard time-advance techniques. 
With the pressure eliminated and continuity satisfied, only two unknowns remain. 

2.3.2 Velocity representation 

Given the weak formulation 2.15, all that remains is to select the basis vectors to 
represent v and £. Using periodicity in 0 and z, allows the following representation 
to be adopted 

Nr 

v(r,0,z,t) = Y^^2Yl a jrni(t)ui(r; k 6y k z )e ,{k6e+kzZ) (2.16) 

j m /-i 

z ) = w/»(r; kg, k z ')e~ i( - ke ' s+k/z ^ (2.17) 

where u / are the basis expansion vectors and W/> the basis weight vectors which are 
both function of r, kg and k z , a jm i are the unknown (complex) expansion coefficients 
and 


27rm 

k z = -j — , — A r e/2 <k e =j< N e /2 - 1, -N z /2 <m< N z / 2 - 1 

k z and kg (or j) are respectively the streamwise and azimuthal wave numbers, and 
L z is the period, or domain length, in the z direction. 

Because v is a real valued quantity, only half of its Fourier coefficients need to be 
computed since there exists a symmetry of the coefficients such that = v* ml . 

From the standpoint of using this symmetry property, the choice of azimuthal wave 
numbers is purely arbitrary; instead, one could easily use symmetry of the axial wave 
numbers. But for the purpose of applying the regularity conditions (see below), using 
symmetry of the azimuthal wave numbers is the natural choice. 

Substituting 2.16-2.17 in 2.15, the discrete form of the equations is obtained 


Act = Bet + f nl 


(2.18) 
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where A and B are respectively the N r x N r mass and viscous matrices, defined as 


rFt 2 

A = { ai'i } = / w/» • u / rdr 

Jo 


(2.19) 


B = {bi'i} = — — / W/' • (V • Vui] rdr (2.20) 

Rei, Jo K ' 

where volume integrals were reduced to single integrals in r by the orthogonality 
property e im4> e~ in4> dd> = 2% S mn . Clearly, the row index corresponds to the weight 

vectors w;», and the column index to the expansion vectors U/. Integrating the viscous 
matrix by parts yields 


B = 


1 

Reb 


/: 


R? — dui 

Vwj- : Vu/ rdr — R 2 \Vr{R 2 ) • ~q^ 


( 2 . 21 ) 


K 2 J 


where a : b is the double contraction operator (= aijbij in tensor notation) and 
V is the Fourier transformed del operator, fni is a vector of length N r accounting 
for the nonlinear term. 2.18 is a linear set of ordinary differential equations for all 
combinations of k$, k z wave numbers. 


Regularity conditions 

Because of the coordinate singularity that appears in the V (= e r + 

JL ej,) operator when written in cylindrical coordinates, special care must be taken to 
ensure that any quantity evaluated at the origin be bounded and smooth. Regularity 
conditions imply a set of requirements imposed on the velocity field; with the velocity 
field being C °° , upon a transformation from cylindrical to cartesian coordinates, the 
origin of the cylindrical system should map smoothly as any other point would. The 
origin is not a special point. 

Those conditions are the subject of appendix B. It is shown that regularity of the 
velocity field amounts to essentially two conditions: first, the radial and azimuthal 
velocity components are constrained to each other. Second, each velocity component 
must have a strictly defined behavior in r near the origin. 
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Vector shape functions 


Within the confine of the regularity conditions, there is a great deal of freedom in 
constructing the expansion vectors. By choosing w // to be the complex conjugate of 
u*, a Galerkin method is obtained in 2.15. Galerkin’s method has the advantages of 
minimizing the error in the L 2 norm and generating a positive definite mass matrix 
(see for example Stanaway et a/.[45]). 


With continuity satisfied and the pressure eliminated, only two degrees of freedom 
are associated with each Fourier/ b-spline mode. Following the notation of Leonard 
and Wray [20], it is convenient to divide the expansion and weight vectors, U; and 
w //, into two distinct classes of vectors (u/ + and U; _ ) and (w*# + and w //“), with co- 
efficients af ml (t) and respectively. The following vectors meet both regularity 

conditions, provided that g x satisfy the appropriate conditions at r — 0 


u i + (r;k e , 



fu r f\ 


( 0 \ 

/ 

k z ) = 


= V x 

o 

= *, 


K u z f / 


s- 

N 

1 

l 



< -*g, > 


( ~ik z g l \ 

u / = V x 

9, 

— 

k z g x 


0 ) 


W+^sJ 


-*hg, \ 

r 9,' + 9 , 

0 / 


( 2 . 22 ) 


(2.23) 


where Vx is the Fourier transformed curl operator, g t = g,(r) are the b-splines 
polynomials, and g/ = dgjdr. The above representation is incomplete when k z = 0, 
i.e. v r — ve = 0; in this case, the k z factor is removed from 2.22 and 2.23 is left 
unchanged. 


Because b-splines are polynomials of any order, they generally do not have the 
correct behavior in r near the origin. To remedy this, b-splines must be made to have 
a behavior analogous to the requirements B.8-B.10. The procedure by which this is 
achieved is the subject of appendix B. A summary of the expansion vectors is given 
in table B.l. 


Upon splitting the vectors in two classes, the velocity becomes Vj m / = aj ml u/ + + 
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a Jtn/ u/ _ and 2.18 is composed of the following 


A+d + + A + d = 
Ala - + A+d + — 

where 


|a + + B + a + f„i + 

(2.24) 

Zct + Btar + + f„f 

(2.25) 

rR 2 

A + = / w i~ • U; + r dr 

(2.26) 

Jo 


f Ri 

A + = j w / ,+ u i r dr 

(2.27) 


r«2 

At = / w// + -U; + rdr 

Jo 

fR-2 _ __ 

AI = / w // • u/“ r dr 

Jo 


and similarly for the viscous matrices. Appendix C gives the complete details on how 
the discrete system of equations is assembled and each matrix is defined. 


2.3.3 Nonlinear Term 

As it appears in 2.24 and 2.25, the nonlinear term is defined as 

ci* = { /.w = S-V- r r r ^ ^ >< r * <» * <2.28) 

This term is computed in the standard way using the so-called pseudo spectral 
approach, where the product w;/ ± ■ (v x «) is computed in physical space by means 
of fast Fourier transforms (FFT) and the final result given by the inverse transform 
(the integral ^77 fo* So*’" e~^ ke ' e+kl '^ dd dz is the very definition of the inverse 
Fourier transform). By using FFT, the pseudo spectral method has the advantage 
of requirering fewer operations than the full spectral approach: typically, the pseudo 
spectral approach requires ~ N log 2 N operations, whereas the full spectral approach 
requires ~ N 2 operations. Fourier transforms are done using the 3/2 rule to avoid 
any aliasing errors. Problems related to aliasing errors have been well documented 
and are not repeated here (see Rogallo [37], Canuto et al. [7] and Sondergaard [42] for 
background). Numerical errors in the radial direction are alleviated by computing all 
integrals to machine accuracy by using Gauss quadrature technique. Specific details 
on how this procedure was implemented are given in appendix C. 



18 


CHAPTER 2. METHODOLOGY 


2.3.4 Time Advance 


The time advancement method used to compute the time derivative in 2.18 is the 
so-called SMR method (see Spalart et al. [44]), which is a mixed explicit-implicit 
method; it uses a third order Runge-Kutta (RK) method for the explicit part and 
a trapezoidal (second order), or Crank-Nicholson, method for the implicit part. An 
equation to be time marched is split into a linear and a nonlinear part 

< 9 $ 

— = L{<$>) + N{*) (2.29) 


where L($) is the linear term which consists of the viscous term and N($) is the 
nonlinear convective term. The method to advance 4>„ at time t to 4> n+1 at time 
t -f At is given by 


<*>' = $ n + A<[L(a 1 $ n + A$0 + 7i^n)] (2.30) 

= <*>' + At[L{a 2 <f>' + f3 2 <t>") + 7 2 N (&) + (! N(<P n )} (2.31) 

$ n+1 = <*>" + Ai[L(a 3 $ ,/ + ^n + i) + 73A r (^") + C2^($ / )] (2-32) 



ai = 96’ 

° 2 40’ 

«3 = — 

6 



B _ 37 
13 160 ’ 

II 

CN 

0 


8 

15’ 

to 

II 

tol 01 

3 

13 = J, Cl 

II 

1 

OD I i — * 
O | -4 

C3 = -^ 


Even though this method requires two previous levels of storage, it can be made 
to require only one previous level of storage (details are in given in section C.4). 


Stability and Accuracy 

The advantage of treating the viscous term with an A-stable implicit method is to 
relieve the time constraint, since it is customary to ignore the implicit part of the 
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equation when establishing that constraint. Consider the following linear model prob- 
lem 


V + u • V<F = 0 

at 


(2.33) 


Assuming that 0 is periodic in all three directions, 2.33 becomes 



u r k T + ug h u z k z ) $ = 0 

k r / 


(2.34) 


For the third order Runge-Kutta scheme the CFL number is \/3, and the stability 
criterion is given as 

At < (2.35) 

| ^max | 

with 


|Amax| = max 
using the triangle inequality 


< max 

< max 


kg 

u T k r + Ug 1- u z k 2 

r 


kg 

|u r |£ r + |u#| h \u z \k 


KIt- + Kl— + Kit 


'A r 


(2.36) 

(2.37) 

(2.38) 


where fc* max = Ng/2 - 1 and & 2max = fi(N x / 2 - 1) and k Tm „ = = f- r with 

A r = R 2 /N r . 

Since maximizing A guarantees stability, assuming periodicity of the radial di- 
rection yields a conservative estimate. Figure 2.2 shows a plot of the modified (or 
numerical) wave number versus the exact value. It is clear that the b-spline approxi- 
mation, in the context of a Galerkin method, yields wave numbers always less than, 
or equal to the exact value. 

Having ensured the solution is well behaved at the origin (through the imposition 
of the regularity conditions), a problem still remains near the origin. This problem 
is best explained by considering a grid in physical space. In the example of figure 
2.3, the mesh on the left shows a concentration of grid point as the radius decreases; 
this would lead to a highly restrictive time step since the value of r A0 becomes very 
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Figure 2.2: Modified wave number ( k ) versus wave number ( k ) normalized by & max , 
for splines of orders 1 through 4. The straight line is the spectral limit (exact) and 
higher order are closer to spectral. 

small near the origin. To alleviate this, grid points should be removed to keep the 
value of r A# approximately constant. Since no special feature (such as a boundary 
layer) exists near the center, there is no need for the added resolution. 

In the present context, i.e. using a Fourier representation for the azimuthal di- 
rection, such a procedure would translate into letting the maximum azimuthal wave 
number vary with the radial position, such that the ratio kg m&x /r appearing in 2.38 
be approximately constant. Because b-splines have local support, splines for larger 
values of kg can be removed near the center of the pipe, while leaving the ones on 
the outer region. This procedure is termed modal reduction. In order for the ratio 
kg max/ r to be constant, kg has to be made a function of r*; the following was adopted 

^max( r ) = (^max “ &)tT + b for 0 < T < R 2 (2.39) 

t 1-2 

t Making k$ a function of r complicates the computation of the nonlinear term since Fourier 
transforms must be done on fixed lengths. Even though this is not a problem in principle the 
implementation of modal reduction in the context of FFT is not trivial. Details are given in section 
C.2 
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Figure 2.3: Mesh (left) with constant number of grid points on each radial location. 
Mesh (right) with constant r AO between grid points. 


where kg max = Ng/2 — 1 and b is set to some small value (2 or 3). At r = 0 there is 
only a single non-zero Fourier mode (kg = 0 for u z and k$ = 1 for u T and ug). The 
modal reduction algorithm can be written as 

otjmi , =0 j > kg m3tX (ri) + 1 (2.40) 


and h are splines with support on r < r,-. Finally, the CFL condition is given as 


At < 


max,- 




x, + l><«, I + Kilfcm.*) 


(2.41) 


where all velocities are evaluated in physical space. The time step can be evaluated 
while the nonlinear term is being computed since it too requires velocities in physical 
space. 


Running with the full CFL number was found to be both stable and accurate. 
In addition to alleviating the CFL condition, modal reduction greatly reduces the 
computations; typically, when modal reduction is active, the computations take 30% 
less time than with the full, unreduced modes. 
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2.4 Computer Code 

The method described in this work was coded to run on the vector-parallel CRAY 
C90 super computer. The computer used had 64 MW § of available memory and eight 
processors. To save memory, the code was designed to take advantage of the solid- 
state storage device (SSD), a very high speed ram disk. The code was written in 
VECTORAL, a computer language designed by Wray [51] to ease implementation of 
vector processing and facilitate the handling of large datasets. Two versions of the 
code were written: one serial (single CPU) and one parallel (some subroutines of the 
parallel version are written in C). Parallelization was straightforward in the context 
of the numerical method and the shared memory architecture of the computer. Both 
are highly optimized to take advantage of vectorization. The serial version runs at 
600 MFlops and the parallel runs at roughly 500 MFlops/CPU, depending on system 
load. Running a grid of 72 x 160 x 192 (N r x Ng x N z ) with fourth order splines, 
the serial code needed 16 MW of memory. For sake of comparison, the code runs 
at 14 jisec/mode/RK sub-step when using cubic splines (without modal reduction), 
whereas the mixing layer code of Spalart et al. [44] achieves 13 /xsec/mode/RK sub- 
step using Jacobi polynomials instead of b-splines. 

The code is divided in two principal subroutines called passes: the first pass 
consist in computing the nonlinear term; it accounts for 85% of the computations. 
The second pass advances the solution to the next time step. The order of execution 
is given as follows 

• Pass 1 

1. Unload the data from database (SSD). 

2. Reorder the data and zero the oddball wave number ( k z — —N z /2) in 
preparation for the fast Fourier transforms (FFT). 

3. Compute Ai~Ki (see section C.2 for the definitions of those terms) and 
transform (FFT) to physical space. 


§MW = mega word, 1 word = 8 bytes. 
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4. If first RK sub-step, compute the maximum eigenvalue and update the 
time step (A t). 

5. Compute 7 t t, $7, Y,7, YV and S t > and FFT to Fourier space. 

6. Obtain the result fnlt and f n i;, . 

7. Save the result to the database. 

8. Proceed to pass 2. 

• Pass 2 

1. Unload the result from pass 1 and the previous step variable from the 
database. 

2. Impose the regularity and boundary conditions on the results of pass 1. 

3. Compute the mass and viscous matrices (A and B). 

4. Compute the effective mass and viscous matrices (A and B). 

5. Impose the regularity and boundary conditions. 

6. If first RK sub-step, do some run-time diagnostics. 

7. Time march (see section C.4 for the exact details). 

8. Save the results to the database. 

9. Proceed to pass 1. 


2.5 Boundary and Initial Conditions 

2.5.1 Boundary Conditions 

The no-slip boundary conditions, where v = 0 on r = ( see 2.8) is used to simulate 
pipe flow. Applying those conditions to the functions given in table B.l, constrains 
the expansion coefficients which have support on r = /?2- The resulting conditions 
are given in table 2.1. 
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Table 2.1: Expansion coefficients for the no-slip boundary conditions. 



k e = 0 

k e >0 

t* 

ll 

o 

<*k — o 

a N r = 0 

a tf r - 1 = a AV = ^ 

a N r - 1 = a A'r = 0 

k z ^ 0 

< = o 

Q Nt- 1 = a Nr = 0 

a N r - 1 ~ Q N T = 0 

Q Nr - 1 " a Nr ~ 0 


In all cases, the support rules given in section A. 3 were used to derive the con- 
ditions of table 2.1. A consequence of those boundary conditions is for = 0, 

which is required by continuity. 


l/?2 


2.5.2 Initial Conditions 

As mentioned in the introduction, pipe flow is linearly stable to infinitesimal distur- 
bances; the onset of turbulence can only be brought about by nonlinear interactions. 
This implies the flow cannot be started with unstable modes in order to trigger tran- 
sition and bring about the onset of turbulence as rapidly as possible. The flow was 
started with a mean velocity profile, or (0,0) mode, given by 1 — r 8 , and with the 
rest of the modes filled with noise with an amplitude of & 10% of the mean. This 
amplitude was found to be more than adequate to excite the nonlinear process. 


2.6 Validation 

Three different approaches were used to validate the code: representation tests using 
known velocity fields, solution to Stokes’ problem and the propagation of linear Orr- 
Sommerfeld waves. Each approach serves a different purpose and is explained below. 

2.6.1 Representation Tests 

The objective of this first test is to ensure that the method was implemented correctly 
(validation of the method proper is left for subsequent tests). The idea behind this 
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test is quite simple: use an arbitrary function instead of b-splines representations 
in the expansion vectors (table B.l) to generate a benchmark result; then 
check that solution against the one obtained by b-splines. The coefficients (o^) 
used for the b-splines solution are generated using the projection procedure shown 
in A.1-A.5. Those tests can be done to verify the mass and viscous matrices, and 
the nonlinear term; they are repeated for all four families of expansion vectors. The 
arbitrary function should satisfy the boundary condition and be regular (possess the 
correct power of r). As an example, one could use 

f(r-k,X) = r'‘«-l(l-r ! ) 2 (2.42) 


where it is assumed R 2 = 1. 


2.6.2 Stokes Flow 

Once the implementation of the method is validated, we turn our attention to the 
validation of the method itself. To validate the method, we make use of Stokes’ prob- 
lem for which there is an exact solution (see appendix D). The objective is to solve 
numerically Stokes eigenvalues and eigenfunctions and compare them to the exact 
solution. The relevance of solving Stokes eigenproblem was best explained by Hughes 
[13] (p. 434): “The accuracy of the eigenvalues and eigenfunctions are measures of 
the quality of both M and K *, in other words, the entire spatial discretization.” 
Hughes [13] gives the following error estimates for eigenvalues and eigenfunctions 
computed using Galerkin’s method when solving an elliptic eigenvalue problem (of 
degree one) 

A s - A, < c/i 2S Af +1 (2.43) 

|| v s -u* || 0 <c/ l s+, A< s+1 >/ 2 (2.44) 

where h is the mesh size parameter, c is a constant independent of h or A 5 , || * || 0 is the 
L 2 norm and 0 < A! < A 2 < . . both A* and v 5 are the numerical approximations. 


^i.e. A and B, the mass and viscous matrices 
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Figure 2.4: Convergence plots for the Stokes eigenvalues (a) and eigenfunctions (b) 
on an equispaced grid with Nk knots, for kg = 2, k z = and s = 5. Plots are for 
cubic x, quartic o and quintic + b-splines. 

Figure 2.4 gives relative errors for the fifth eigenvalue and eigenfunction for kg = 2 
and k z = 27 t ; plots for other combination of wave numbers are included for complete- 
ness and given on figure 2.7, page 30. As expected, errors for the eigenvalues are much 
lower then their eigenfunctions equivalent since as shown in 2.43 and 2.44, eigenvalues 
converge faster than eigenfunctions (exponent of h ). Eigenvalue errors quickly reach 
round-off levels for higher degree splines and moderate values of Nk (the number of 
knots). 

Slopes of each curve represent the convergence rates and are given in table 2.2. 
Computed rates are always bounded from below by the estimates 2.43-2.44, and the 
agreement between them is very good. Preserving the full rate of convergence gives 





2.6. VALIDATION 


27 



Figure 2.5: Relative error in eigenvalue for N T = 30, k$ = 0 and k z — 2it. Symbols 
are as figure 2.4 

a strong indication that the computations are behaving as expected. 


Table 2.2: Convergence rates for k$ = 2, k z = 2ir. 



Eigenvalues 

Eigenfunctions 

5 

computations 

(2.43) 

computations 

(2.44) 

3 

4.14 

4 

3.11 

3 

4 

6.50 

6 

4.21 

4 

5 

9.34 

8 

5.42 

5 


Because the expansion functions as defined in table B.l involve derivatives of b- 
splines, with a b-spline of degree 5, the method will converge at the rate of a spline of 
degree 5 — 1; this should not be surprising since the g th derivative of a b-spline a degree 
5 can be represented by a b-spline of degree S — q. This explains the discrepancy one 
might observe between the estimates 2.43-2.44 and the values of table 2.2. 

This behavior is best observed in the expansion functions for k$ = 0 and k z > 0: 
the u/ - vector involves a derivative of a b-spline whereas U/ + does not. For this case, 
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the plus and minus modes should converge at different rates. Figure 2.5 clearly shows 
this behavior: the plus mode of the cubic spline converges at the same rate as the 
minus mode of the quartic. Figures 2.7 (e) and (f) show that for the (0,0) modes, 
convergence is much faster since these expansion functions involve no derivative of 
b-splines. So, every derivative incurs a loss of an order of accuracy. Figure 2.5 also 
shows an interesting consequence of the estimates 2.43-2.44. Because the eigenvalue 
A s appears on the right hand side of the estimates, when the eigenvalue becomes large 
the error follows accordingly; with s large, the error becomes 0(1). 

2.6.3 Wave propagation 

The last test consists in propagating Orr-Sommerfeld waves. The objectives are to 
validate the time marching method, and the computation of the nonlinear term. It 
also ensures that the code as a whole is working properly, such as access to the SSD 
and other I/O. 

The idea is to start the computations from a solution to the linearized Navier- 
Stokes equations for an arbitrary wave number pair (£#,£*); such solution is known 
as an Orr-Sommerfeld wave (see appendix D) and is given by D.19. This wave is then 
propagated using the standard Navier-Stokes solver. If the amplitude of this wave is 
small enough, it should behave linearly such that the solution at any given time is 

u (r,t;kg,k z ) = eu(r,0-,ke,k z )e~ %u,i (2.45) 


where e is a small number chosen to be 10 6 . 

Figure 2.6 gives the error norms for the four different wave number pairs after 
having propagated for fifty time steps. The linear growth of the error in time for 
plots (a)-(c) is in agreement with the predicted error, which has a form 

\e(t n )\ = cAt 2 t n + 0(At 3 ) (2.46) 

where c is a constant independent of At, and t n — nAt. 2.46 implies that if At is 
small enough, so higher order terms become negligible, the error will grow linearly in 
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Figure 2.6: Relative L 2 error norm of the Orr-Sommerfeld wave for the four different 
wave number pairs (k$,k z ): (a) (1,1); (b) (1,0); (c) (0,1) and (d) (0,0). Parameters 
are Re b = 9600, N r = 35 ,S = 3, N 8 = N z = 12 for a CFL = 0.05. 


a finite time. In absolute terms, the low error is also an indication that the nonlinear 
term is behaving as expected by preserving linearity of the solution. Plot (d) is a 
simple measure of how well the mean flow is preserved, since the nonlinear term has 
zero contribution for the (0,0) mode; the error is no larger than round-off. 

Lastly, it was possible to observe an additional wave or alias at (2k 8 ,2k z ) with 
amplitude 0(e 2 ) being generated by the nonlinear term. This is not unexpected since 
the product v x w will generate a term with an exponential of the form e 2fc <?+ 2fc *. This 
is further evidence that the nonlinear term is behaving as it should. To make this 
alias disappear, we could have selected a smaller initial amplitude such that e 2 would 
be of the same order as round-off. 
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Chapter 3 
Pipe Flow 


We present in this chapter results for the pipe flow simulation. The following results 
were obtained using the method described in chapter 2. The pipe flow simulation 
of Eggels et al. [12] (hereinafter referred to as EUW), the channel flow simulation 
of Kim, Moin and Moser [17] (KMM), the PIV (Particle Image Velocimetry) and 
LDA (Laser Doppler Anemometry) measurements in pipe flow of Westerweel et al. 
[50] (WDV), and the high resolution LDA measurements (also in pipe flow) of Durst, 
Jovanovic and Sender [11] (DJS); all were used to compare with the present results. 


3.1 Test Case 

In wanting to keep within the same Reynolds number range as the aforementioned 
investigators, the Reynolds number based on the bulk velocity, Ub, and pipe diameter, 
D, was set to Reb = UbD/v = 5600 (or a Reynolds number of Re T — 380 based on 
the wall shear velocity u T ). This value matches KMM’s value of 5600 based on the 
mean velocity and channel width (or Re r = 360), and is close to EUW’s value of 5300 
(Re T = 360) based on Ub and D. The computations were carried out on a grid of 
72 x 160 x 192 ( N r xJV«x N z ), for a total of 2.2 million Fourier/b-splines modes (less 
if we account for modal reduction) with quartic b-splines. The length of the domain, 
L z , was chosen to be 5D, the same as EUW; the adequacy of the spatial resolution 
and L z will be assessed later by examining the velocity and vorticity spectra, and the 
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streamwise two-point correlations. With this computational domain, grid spacings 
are about /?2A0 + ~ 7.5 and Az+ ~ 9.9 wall units*. A non-uniform grid is used in 
the radial direction based on an exponential function: close to the wall the radial 
grid spacing is smaller, with the first point away from the wall located at r + rs 0.39; 
spacing reaches its maximum near the center with Ar + rs 5.7. Very close to the 
center, a finer radial grid is used in order to impose the regularity conditions; a finer 
grid makes this process more local. 

Those values compare well with EUW who used 96 x 128 x 256 grid points with 
an equispaced mesh in the radial direction. The first point away from the wall was 
located at r + = 0.94, with the rest spaced by Ar + ~ 1.88. Other grid spacings were 
-/?2A0 + w 8.84 (at the wall) and A z + ~ 7. Since EUW do not remove azimuthal 
modes as the radius decreases, the azimuthal resolution varies linearly with r and 
reaches a minimum of 0.05 wall units near the centerline. Even though the total 
number of grid points exceeds the one used for our present simulation, EUW’s res- 
olution must be considered at par if not inferior to the present one: since EUW’s 
computational procedure is based on second order finite differences, which possesses 
a high diffusivity compared with the present technique, many more grid points must 
be used to adequately represent the small scales. EUW’s choice of a uniform grid in 
the radial direction is certainly questionable. 

KMM’s computations used 192 x 129 x 160 modes in the streamwise, crosstream 
and spanwise directions (y,y,z) respectively. KMM’s grid spacings were Ax + ~ 12 
and A z + rs 7. The channel had a length of 2nD (versus 5 D for the pipe) and a 
spanwise length of ttD, where D is the channel width, or distance between walls. The 
crosstream, or y direction, was discretized with a cosine grid (necessary for Tcheby- 
shev approximation) with the first point away from the wall located at y + « 0.05, 
reaching a maximum spacing of 4.4 at the center. Their near wall resolution is almost 
an order of magnitude smaller than the value adopted for the present simulation. 
There appears to be very little advantage in using such a fine resolution since the 


"the -f superscript refers to nondimensionalization with respect to wall variables, mainly the wall 
shear velocity u \ = (t w / p) = — u du z /dr e.g. r + — ru 7 /v. 
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viscous sublayer is usually of order y + = 5^. The only apparent reason for this, is an 
attempt to resolve numerical artifacts produced near the boundary (so-called projec- 
tion error ; see Moser and Moin for background [25]). Because of the cost involved 
and the lack of additional physical insight gained by such a fine grid, no attempt was 
made to duplicate their resolution. 


3.2 Turbulence Statistics 

The results presented in this section are comprised of 46 different fields approxi- 
mately equispaced in time, averaged over a period of 43 time units (based on the 
non-dimensional time D /Ub )• The results of EUW consist of 41 fields averaged over 
a period of 59 time units. Although 40 time units might not be ideal to remove all 
statistical errors, it is certainly adequate. Statistical steady-state is reached when 
the total shear stress t tz — 1 /Rebdu z /dr - u' T u' z becomes a linear function of r. If 
r rz is normalized by t w , then dr+fdr - 1. Defining the deviation from linearity as 
|| T r + — r||o/||r||o, the deviation for the present simulation is 0.51%, and for the DNS 
of EUW, 0.72%. 

3.2.1 Mean flow 

Table 3.1 compares several mean flow parameters reported by various investigators 
to the results from the present computation. Of interest, is the friction coefficient 
defined as 

C, = 2(g) (3,) 

The computed value of 9.16 x 10 -3 agrees within half a percent of Blasius’ empirical 
formula. Differences between the present simulation and EUW values for Cj , the dis- 
placement thickness 8 , and the momentum thickness 0, are consistent with expected 
dependence on Reynolds number (with the Reynolds number increasing, those val- 
ues should decrease). Differences in geometry between the pipe and channel are also 


tFor the pipe, y + is defined as y + = (R 2 — r)u T /v. 
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Table 3.1: Mean flow properties from several investigators. Repeated in modified 
form from EUW. 



b-splines 

DNS (EUW) 

PTV (WDV) 

LDA (WDV) 

HWA (EUW) 

DNS (KMM) 

Re c a 

7248 

6950 




6600 

Reb 

5600 

5300 

5450 

5450 

5600 

5600 

Re T 

380 

360 

366 

371 

379 

360 

U C /U T 


19.31 


19.39 

19.40 

18.20 

Ut/u T 


14.73 

n 

14.68 

14.76 

15.63 

Uc/Ut 

1.29 

1.31 

1.30 

1.32 

1.31 

1.16 

C f (x 10" 3 ) 

9.16 

9.22 

9.03 

9.28 

9.18 

8.18 

wammmimi 

9.13 

9.26 

9.19 

9.13 

8.44 c 

EHBEBUI 

9.36 

9.51 

9.43 

9.36 

- 

6 /R 2 ' 

0.121 

0.127 

0.124 

0.130 

0.128 

0.141 

0 /Rj 

0.066 

0.068 

0.068 

0.071 

0.070 

0.087 

II 

^t> 


1.86 

1.83 

1.83 

1.82 

1.62 


webm 

8.91 

8.78 

8.79 

8.73 

6.97 


“All Reynolds numbers are based on D, the pipe diameter or channel width, and v. 
J Friction coefficient based on Blasius’ formula C/ bU = 0.079Re^° 25 . 
c Based on Dean’s correlation Cj = 0.073 Re^° 25 (from KMM). 

^Friction coefficient based on Barenblatt’s hypothesis (see page 38). 

*The displacement thickness 6 is defined as 6(2i?2 — <5) = 2 0 ~ d z (r)/U c ) rdr (EUW). 

^The momentum thickness 0 is defined as 0(2i?2 — 0) = 2 f 0 2 u z (r)/U c (\ — u z (r)/U c ) rdr. 
9G = U c /u T {{H -l)/H). 


apparent; for a similar Re T , the channel has a significantly lower friction coefficient, 
with larger 8 and 6. 

Figure 3.1 shows the mean velocity and vorticity profiles. Agreement with EUW 
and the experiments of WDV is excellent. The slight bulge in the velocity profile for 
0.3 < r < 0.8 compared to EUW is due to a slightly higher Reynolds number in our 
case, which causes the velocity profile to be fuller. As expected, the mean vorticity 
is maximum at the wall and drops to zero at the center, indicative of a vortex sheet 
in the near-wall region. 

Figure 3.2 (a,b) shows the mean velocity profile normalized by the wall shear stress 
velocity compared to the experimental data of WDV, and the results from EUW and 
KMM simulations. Figure 3.2 (a) reveals good agreement between the experiments 
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(a) 



(b) 



Figure 3.1: Mean profiles: (a) stream wise velocity normalized by the centerline ve- 
locity; (b) azimuthal vorticity normalized by the vorticity at the wall. 



Figure 3.2: Mean velocity profiles scaled with 
the pipe and channel. 
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and the two simulations. In the viscous sublayer ( y + < 5) the present simulation 
and the law of the wall collapse to a single curve until y + ~ 4. This collapse is 
even more apparent in figure 3.2 (b) where the pipe is compared to the channel; all 
three curves are undistinguishable in the viscous sublayer. In the inertial sublayer 
(log region), the expected discrepancy between pipe flow and the log law is apparent. 
Patel and Head [32] observed this discrepancy for pipe flow up to Ret, ~ 10000 + . A 
closer look at figure 3.2 (a) shows that in the inertial sublayer, the present simulation 
is closer to the logarithmic profile than that of EUW, which is consistent with the 
trend observed by Patel and Head for higher Reynolds numbers. For channel flow 
however, they found that a lower value of iZej « 3000 was sufficient to match the 
logarithmic profile. Figure 3.2 (b) compares the velocity profiles for the pipe and 
channel; following KMM, the additive constant is taken to be 5.5 to compensate 
for low Reynolds number effect, instead of the more common 5.0 (see Tennekes and 
Lumley [47] chap. 5). As expected by Patel and Head’s [32] prediction, the channel 
matches the logarithmic region. 

Recently, Barenblatt [2] and Barenblatt and Prostokishin [3] introduced a power 
law to describe the velocity profile in the inertial sublayer instead of the traditional 
logarithmic profile. The novelty in Barenblatt’s approach is the dependency of the 
predicted profile on the Reynolds number. Whereas the log law is independent of the 
Reynolds number (based on Karman’s constant k — 0.40 and an additive constant of 
5.0), Barenblatt proposed the following 

8+ = (4= In Re + s ,+ 3 « 21 " s ') (3.2) 

where the Reynolds number can either be based on the bulk or centerline velocity 
(this was found to make little difference even in the low Reynolds number range). 
Figure 3.3 compares the velocity profiles for the log law and Barenblatt’s power law 
(we should point out that Barenblatt makes no claim as to the validity of his approach 


1 Although the recent experiments on pipe flow by Durst, Jovanovic and Sender [11] at Ret, = 7442 
revealed perfect match with the logarithmic profile, raising doubts as to what the exact threshold 
should be. 
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Figure 3.3: Comparison of velocity profiles in the inertial sublayer with the standard 
log law and Barenblatt’s power law. 


in the low Reynolds number range). The power law underpredicts the velocity in most 
of the range of y + . Because Barenblatt’s power law is Reynolds number dependent, 
it was possible to suggest a skin friction law given by 


where 


C 


I ba 


2 

^2/(l+o) 


e 3/2 (v/3 + 5a) 

2 a a(l + a)(2 + a) 


and 


3 

2 In Re 


(3.3) 


(3.4) 


Results are given in table 3.1. It is clear that the agreement is far from opti- 
mal, with errors in the 2 to 3% range, although the error seems to diminish as the 
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Reynolds number increases. This would be consistent with Barenblatt and Prostok- 
ishin’s [3] observations that for higher Reynolds numbers, the predicted values are 
within experimental scatter (although they only used Nikuradse’s rather dubious pipe 
flow measurements; from figure 3.3 Nikuradse’s measurements recover the logarith- 
mic profile even though for this Reynolds number range, the actual profile should lie 
above the log region). In summary, Barenblatt’s power law approach does not pro- 
duce better results than the standard logarithmic profile in the low Reynolds number 
range, although the method yields close estimates of friction coefficients. 


3.2.2 Spectra and two-point correlations 

Assessing the resolution of the computational domain can be made through close 
examination of the one-dimensional velocity and vorticity spectra. Figures 3.4-3.5 
and 3.6-3. 7 respectively show the velocity and vorticity spectra for several radial 
locations: one in the viscous sublayer (j/ + = 3.458), one in the buffer layer ( y + = 

21.99), one in the inertial, or log region ( y + = 122) and one the wake region ( y + = 

169.1). The one-dimensional velocity spectra are defined as 

£«(r; k.) = 4- J r Q«(r, sey^'dff (3.5) 

£«( r; k,) = C Q (c (r,6z)e ii "' iz' (3.6) 

L z Jo 

where 

Q^(r, SO) = u' ( (r, 0, z, t)u' ( (r , 0 + 6 ' , z, t ) (3.7) 

Q a {r, Sz) = «£(r, 0 , 2 , t)u[{r, 0,z + z', t ) (3.8) 

are the two-point correlations for streamwise and azimuthal separations (no sum on 
repeated Greek indices). The velocity spectra is simply the Fourier transform of the 
two-point correlations and vice versa. Vorticity spectra (0^) are similar to 3.5 and 
3.6, with u>', replacing u £. 

Velocity spectra reveal that all small scales are adequately resolved since no sig- 
nificant accumulation of energy is observed at high wave numbers. Away from the 
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Figure 3.9: Comparison of normalized velocity spectra with LDA results for y + ~ 11. 
B-splines , LDA *. 


wall, smaller scales show at least four decades of reduction of energy compared to the 
larger scales (low wave numbers). Near the wall, significant anisotropy is detected 
with energy of the radial perturbations being almost three orders of magnitude less 
than the streamwise perturbations at lower wave numbers; this is especially apparent 
in the azimuthal spectra. Conversely, the wall has a sustaining effect on the stream- 
wise perturbations by maintaining a higher energy level at low wave numbers. Large 
radial perturbations are virtually non-existent near the wall. Figure 3.4 shows energy 
spectra compared to the inertial range of the Kolmogorov spectrum, or k 5//3 law. 
Away from the wall (typically in the log region), The energy follows the Kolmogorov 
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spectrum rather closely at moderate wave numbers, indicating mostly turbulent dif- 
fusion; closer to the wall viscosity becomes predominant and the energy decays faster. 
The higher level of noise in the spectra near the center are caused by the small number 
of azimuthal modes, which degrades spatial averaging. 

Further evidence of the adequate resolution is seen in the vorticity spectra, since 
vorticity is more sensitive to the resolution of the small scales. Vorticity spectra 
show at least three decades of decay and no significant pile-up of energy at high 
wavenumbers. A very small amount of pile-up is observed in the streamwise spectra 
for y + ~ 22, but overall should have no consequence because of the large drop in 
energy between low and high wave numbers. Near the wall, the higher energy level 
of the azimuthal perturbations is a direct consequence of the streamwise velocity 
perturbations also being large (see figure 3.16 and 3.18). The rather abrupt cut-off 
observed at high wave numbers near the center is a consequence of modal reduction. 

Figures 3.8 and 3.9 compare the present spectra with those from the experiments 
of WDV and the DNS of EUW. The spectra were normalized as follows 


E+(k,D) 


£«(*,£) 

u\k (£«( 0)/2 + Eli?" 1 £«(»*)) 


(3.9) 


where k = 2-kD/ L z . Figure 3.8 shows almost perfect agreement between the exper- 
iments and the present results for k z D up to « 40, after which the PIV results are 
unreliable. Agreement between the two simulations is also good, but two important 
differences appear: first, as expected, the higher wave numbers of the DNS of EUW 
diffuse much more rapidly than the present results; this is consistent with the higher 
diffusivity of their numerical scheme. Second, the near wall streamwise spectra of 
EUW undershoots the one obtained by the b-splines over most of the spectrum (see 
figure 3.8 for y + = 18); this is probably due to the lack of near wall resolution com- 
bined with the higher diffusivity of their method. This shows that most of the smaller 
scales are not adequately resolved by EUW. 

Figure 3.9 is similar to figure 3.8, but shows the higher resolution of LDA tech- 
nique. The differences between the present simulation and the LDS measurements 
are difficult to explain, in light of the almost perfect agreement of the PIV results. A 
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possible explanation is the proximity of the wall ( y + « 11) which makes the measure- 
ments more difficult; an anomaly of the experimental data is clearly visible in figure 
(b) at k z D £« 45. 

Figures 3.10 and 3.12 give the two-point velocity correlations for stream wise and 
azimuthal separations respectively; figures 3.11 and 3.13 give the equivalent two-point 
vorticity correlations. Figure 3.10 for y + ~ 3.5, shows the velocity perturbations are 
still slightly correlated for large streamwise separations which implies that the size 
of computational domain ( L z ) might be marginal (too short), although not by much 
since the values are quite small. For other radial positions, the zero correlation point 
is still within statistical sampling error as indicated by the standard deviation; near 
the center, the error becomes very large due to the fewer number of azimuthal modes. 
EUW’s results also show symptoms of a marginal domain length. However, their 
results show that the correlation is maintained further away from the wall than in 
the present computations; this is again due to their larger viscosity and numerical 
damping. 

The streamwise component of the streamwise velocity correlations are correlated 
over longer distances, indicating the presence of coherent streamwise features, or so- 
called streaks. The two-point vorticity correlation (fig. 3.11) also indicates that in the 
near wall region, the values are correlated over the length of the domain. Note however 
that near the center, the vorticity correlation gives a much clearer picture than the 
velocity correlation did: there is no evidence of correlated values near the center, and 
even down to y + 22. In the viscous sublayer, the presence of long vortex-sheets are 
apparent by the longer correlation length of the azimuthal component. 

The extent of the two-point velocity correlations for azimuthal separations show 
the circumferential size, or width of the streaks. The results show that the width 
of the streaks is rather narrow as indicated by the zero-crossings of the streamwise 
component at y + « 3.5; we obtain A + « 55, where A + is the average width of a 
streak, for an average streak spacing of 2A+ « 110; this is in agreement with Kline et 
al. [18] who experimentally found a mean streak spacing of 2A + « 100 in boundary 
layers, Moser and Moin [25] who also obtained 2A + « 100 in curved channel flow, and 
similarly for KMM and the straight channel. The large negative correlation in the 





Figure 3.10: Two-point velocity correlations for streamwise separation: Q rr , 

Qee , Qzz • The error bars indicate the standard deviation (±<r) of the 

streamwise component. 
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Figure 3.11: Two-point vorticity correlations for streamwise separation. Symbols as 
for figure 3.10. 
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Figure 3.12: Two-point velocity correlations for azimuthal separation; symbols are as 
for figure 3.10. 

radial component is produced by the impingement and ejection of fluid toward and 
away from the wall; this was also observed in KMM and Moser and Moin [25]. Two- 
point vorticity correlation for azimuthal separations are also shown on figure 3.13. A 
large negative correlation is observed which indicates the presence of counter-rotating 
radial vortices. Figure 3.14 shows contour plots of the stream wise and radial vorticity 
and velocity perturbations for a single flow realization; it illustrates physically the 
two- point correlations: the presence of a (short) streamwise vortex (3. 14. a) entrains 
fluid, with one side going toward the wall (impingement, u' r > 0) and one side going 
away from the wall (ejection, it' < 0) (3.14.b). The impinging fluid has a positive 
radial vorticity and the ejecting fluid, negative radial vorticity (3.14. c); the motion of 
the impiging and ejecting fluid is not strictly linear but helical. These two counter- 
rotating radial vortices entrain fluid in the streamwise direction in the region between 
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Figure 3.13: Two-point vorticity correlations for azimuthal separation; symbols are 
as for figure 3.10. 


their respective cores (3.14.d), yielding a high speed streak. The positive streamwise 
velocity perturbation is seen to extend behind the location of the streamwise eddy 
indicating that the eddy is moving downstream in the pipe leaving behind it a high 
speed streak. 


3.2.3 Turbulent intensities 

Figures 3.15 (a) shows the Reynolds shear stress for the present simulation, the exper- 
iments of WDV and computations of EUW. The measurements, especially PIV, show 
significant deviation from the computations; PIV is better suited to collect instan- 
taneous fields than to gather statistics (see EUW). LDA results track the computed 
value more closely, but lack resolution for r > 0.7. Both computations show close 



e and radial velocity and vorticity perturbations ( y + « 5). Length of 
lines are respectively positive, zero and negative contours. 
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agreement, except for r > 0.8 where the differences are consistent with differences in 
Reynolds number. The total shear stress (see page 33), is almost linear indicating 
statistical steady state or fully developed turbulence. The exact value of r+ = r is is 
plotted for reference. The present results agree well except for a slight overshoot for 
0.5 < r < 0.65, whereas EUW has a slight overshoot for r < 0.65. Figure 3.15 (b) 
reveals almost perfect agreement between the pipe and channel. The slight overshoot 
of the b-splines between 0.2 < r < 0.6 is probably caused by the number of samples 
used to generate the average. 

Figures 3.16 and 3.17 show the rms (root mean square) of the three velocity pertur- 
bations. Figure 3.16 (a) reveals relatively good agreement with the experiments and 
both simulations. For the radial and azimuthal components, the differences between 
both computations are consistent with different Reynolds numbers, with the present 
simulation being somewhat higher. The difference for the streamwise component is 
rather surprising, with EUW simulation overshooting the present computation in the 
region 0.6 < r < 0.9, a trend opposite to what is observed for the other two compo- 
nents. The higher energy of the streamwise component is probably due to numerical 
artifacts produced by EUW method, such as inadequate resolution in the near-wall 
region. This is confirmed by figure 3.16 (b) where the same rms velocity fluctuations 
are compared with KMM’s channel simulation. Slight differences between all three 
components are observed, with the pipe results being consistently above the chan- 
nel’s; values are closely matched near the wall, and somewhat higher intensity for the 
pipe near the center, a trend similar to the mean flow. 

DJS conducted a series of LDA measurements in the near-wall region of pipe flow 
at different Reynolds number ranging from 7442 to 20 500. The results for Reb = 7442 
are plotted in figure 3.17; the agreement between their measurements and the present 
computations is excellent, especially for the streamwise and azimuthal components. 
The added noise in the azimuthal component is not surprising since this is the most 
delicate direction to measure (WDV did not even attempt to make measurements in 
the azimuthal direction). In fact, most of the data showing the greatest deviation was 
measured from the centerline to the upper wall, whereas the data measured from the 
centerline to the lower wall agrees much more closely. Measurements for the radial 
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(b) 



Figure 3.15: Reynolds shear stress: (a) pipe data with total shear stress (r r + = r); 
the exact value of r r + is also plotted for comparison, (b) comparison between the pipe 
and channel. 
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(a) 



(b) 



Figure 3.16: RMS (root mean square) of velocity fluctuations: (a) pipe, (b) compar- 
ison between the pipe and channel. 






Figure 3.17: RMS of velocity fluctuations: comparison between the present simulation 
and DJS’ near-wall LDA data (Reb = 7442). 
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component tend to fall between the present results and EUW s simulation; a possible 
explanation for this lower than expected amplitude might stem from DJS’ rejection 
of high amplitude radial perturbations (see section 3.2.4). Farther from the wall, 
Reynolds number effects are apparent with the measurements having higher values 
(figure 16 of DJS shows this trend for different Reynolds numbers). 

Figure 3.18 (a) shows the rms pressure^ fluctuations normalized by pu 2 T . The shape 
of the profile is in good agreement with the results of EUW; the higher values for 
the present simulation is a result of higher Reynolds number, which seems to support 
KMM’s assumption that pressure fluctuations are dependent on Reynolds number, 
even when normalized by wall variables (see KMM section 4.2). The results are also 
in qualitative agreement with KMM, where they obtained a maximum of 1.75 at 
y + « 30 with a wall value of « 1.5, versus a maximum greater than 2 at y + « 31 
and a wall value of ~ 1.65 for the present results. The centerline value is also higher 
at « 0.95 for the pipe versus « 0.75 for the channel. Clearly, for a similar Reynolds 
number the pipe possesses larger pressure fluctuations than the channel. 

Figure 3.18 (b) shows the rms of vorticity fluctuations normalized by u\jv. These 
results are in good qualitative agreement with KMM, Moser and Moin [25] and the 
pipe flow simulation of Zhang et al. [55], the latter being a much lower Reynolds 
number flow (Re b = 2500). Unfortunately, no experimental measurements could 
be found for comparison. Comparing with KMM, the streamwise vorticity shows a 
local minimum at y + ~ 5, with a maximum being reached at the wall. KMM have 
attributed this to the presence of streamwise vortices close to the wall. Figure 3.19 
shows a view in the polar plane of the streamwise vortex shown on figure 3.14. The 
maximum in rms vorticity at the wall is explained by the presence of a mirror vortex 
with vorticity inverse to that of the eddy. Between the mirror vortex and the eddy, 
vorticity reaches zero; in the average sense however, this translates into vorticity 
reaching a minimum. 

For r < 0.75 all three vorticity fluctuations are virtually identical contrary to 
velocity fluctuations which are quite different. This points to isotropy of the small 


§ Because pressure drops out of the equation set (see 2.3.1), a Poisson solver had to be written to 
compute the pressure which can be determined solely from the velocity. 



Figure 3.18: (a) RMS of pressure fluctuatioi 


fluctuations normalized by ul \jv. 
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Figure 3.19: Streamwise vorticity perturbation u)' z . Solid, dashed and dotted lines are 
respectively positive, zero and negative contours. 

scales since vorticity fluctuations are by nature small scales phenomena whereas ve- 
locity fluctuations tend to be produced by larger, more energetic scales. 


3.2.4 Higher order statistics: Skewness and Flatness 


Skewness and flatness (also referred as kurtosis ) are shown in figures 3.20-3.23. They 
are respectively defined as third and fourth moments of velocity perturbations 


Q< >\ U < 

SM = 

U- 


1 3 


- 3/2 ' 


f K) = 


U ’: 1 


u. 


/2 


(3.10) 


Each is related to the probability density: if the skewness is positive, then events with 
large negative values of uf are not as frequent as events with large positive values 
of if the flatness is large, then isolated events with large velocity perturbations 
are possible (see Tennekes and Lumley [47]). For a Gaussian distribution of the 
skewness and flatness are respectively 0 and 3. 

Figure 3.20 compares the skewness obtained from the present computation, with 
the experiments of WDV and the DNS of EUW. Results from the present simulation 
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and PIV measurements show significant disagreement throughout the flow field, at- 
tributable to the low precision of the PIV method. LDA results for the streamwise 
component are in good agreement with both computations except near the wall where 
the results are obscured by noise. The measurements of DJS (figure 3.21) shows good 
agreement except for a slight underpredicting of the the streamwise component for 
y + < 5. In such a close proximity to the wall, accurate experimental measurements 
are hard to achieve. Note that because of the normalization by u( n , skewness and 
flatness results for y + = 0 are excluded, since the numerical limit as the wall is 
approached is undefined, even though an analytical limit might exist. 

Clearly, the most interesting difference is observed between both simulations for 
the radial component. The present results reveal two zero-crossings (r ~ 0.8 and 
r « 0.97), whereas EUW obtained only one zero-crossing at r w 0.75. On the other 
hand, the present results are in agreement with the results of KMM (not shown) 
who also obtained two crossings at approximately the same locations. This seems 
to be confirmed by the LDA results of DJS (figure 3.21) where the trend in the 
data points to another crossing at y + Pa 3, although not decisively. The discrepancy 
between the present results and EUW can only be caused by the coarser near-wall 
resolution of EUW’s computation. Scatter in the azimuthal component is attributable 
to a poor statistical sample, since this component should be strictly zero by the 
equal probability of positive and negative u' e . Similar scatter can also be seen in the 
streamwise component closer to the center. 

A similar comparison of the flatness (figures 3.22 and 3.23) reveals good agreement 
between the LDA measurements of WDV and DJS, the simulation of EUW and the 
present simulation, except in the near-wall region. Typically, for y + < 5 the experi- 
mental measurements become obscured by noise and become unreliable. Coarseness 
of EUW’s grid probably explains the differences between both computations in the 
near-wall region. For example, EUW computed a flatness of 5.5 at the wall for the 
streamwise component, whereas the present simulation gives 4.5; this value compares 
favorably with KMM who obtained 4.2. 

The most important discrepancy is observed between the present simulation and 

















Figure 3.21: Skewness of velocity fluctuations: 
surements (DJS). 
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Figure 3.22: Flatness of velocity fluctuations. F (u() = 3 indicates Gaussian distribu- 
tion. 




Figure 3.23: Flatness of velocity fluctuatio 


surements (DJS). 
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DJS’ measurements of the radial component. Whereas the computed flatness in- 
creases monotonically reaching a maximum at the wall, DJS’ measurements shows 
a maximum at y + ~ 12 and then decreases towards a near-Gaussian behavior at 
the wall. Comparing with KMM’s data, which shows the same trend as the present 
results, DJS speculated that the discrepancy might be caused by the near-wall reso- 
lution of KMM’s simulation, even referring to KMM’s results as “strange”. Another 
explanation for the difference between the simulation and the measurements was given 
by Westerweel (Private communication). He observed that large fluctuations could be 
rejected by the filtering in the LDA processor. Xu et al. [54] conducted experiments 
in pipe flow for the sole purpose of determining whether or not large flatness was 
physical in nature. They were able to determine that events with large radial pertur- 
bations are indeed present in the near-wall region (y + < 5) by taking measurements 
at very short intervals (1 kHz); they also showed that with DJS’ slower sampling (100 
Hz), events with large perturbations could be ignored (thus the filtering effect); in 
addition, DJS interpreted as error and removed any measurement with peak larger 
than seven times the rms velocity. Those large perturbations, or spikes, can exceed 
\u' r ju T rms| = 10, with rather short time scales (20 viscous time units based on u\jv ), 
which make the observations difficult. The present simulation tends to support Xu 
et al. observations: the large radial velocity perturbation shown in figure 3.14 (b) is 
equivalent to a ratio |u(./u r rms| — 7-71 at j/ + = 5. Closer to the wall (j/ + = 0.385, the 
first grid point away from the wall) a ratio of 10.60 is obtained; this also points to the 
relative sensitivity to grid spacing. Recall that with a coarser grid, EUW obtained a 
smaller flatness (see figure 3.22). 


Figure 3.24 (see also figure 8 of Xu et al.) shows the effect of filtering on flatness. 
By removing values larger than a certain factor of the corresponding rms velocity, the 
flatness can be reduced to match the data of DJS. By dropping all perturbations larger 
in magnitude than five times the rms velocity, the results of DJS can be reproduced. 
Note however that DJS use a filter seven times the rms perturbation, but a filter 
of five is found to give the best agreement. This is probably due to the additional 
filtering originating from the slower sampling in DJS’ measurements. 
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Figure 3.24: Flatness of the radial perturbation with different filters: no filter, 

|u r /u rrms | ^5, — . ~ _ |u r /u rrms | ^ 7 . 

Xu et al. attributed the success of DNS in calculating accurate flatness to the na- 
ture of the averaging process: for DNS, at each radial location averages are computed 
on cylindrical surfaces, such as the one shown in figure 3.14, whereas for experiments 
they are measured on lines, which greatly reduces statistical accuracy. 


3.2.5 Budgets of turbulence transport equations 


The transport equations for the Reynolds stresses were derived for cylindrical coor- 
dinates in Moser and Moin [25] ; in tensor form they are given by (e.g. see Mansour 
et al. [22, 23]) 

Jfyi U i U j ~ *b Tij + D;j + II, 'j (3.11) 

where D/Dt = djdt + Ukd/dx^, Uk is the mean velocity and the terms on the right 
hand side of 3.11 are given by 


Pij = - ( u i u kUj,k + u j u kUi,k) 
= 2i/u U u U 

T a = ~ {< u 'i u 'k) k 
D< > = V 

n u = - + u'jP'.i) Ip 


Production 
Dissipation 
Turbulent transport 
Viscous diffusion 
Velocity pressure-gradient 
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where indices after the comma refer to covariant derivatives. A transport equation 
for the turbulent kinetic energy, defined as k = \u\u\ = \{u' r u' T + u' 6 u' 9 + u'u' ) (sum 
on repeated roman indices), can be obtained by taking the trace of 3.11. It is also 
customary to split the velocity pressure-gradient term into a pressure-strain (IISY;) 
and a pressure-diffusion term (IIZ?y), where 

US tJ = 2^- (3.12) 

Sij = l/2(u< ■ + «',) and UD^ = Ily — IlSy. In this form, the pressure-strain term 
redistributes energy among the different components. 

Figure 3.25 (a)-(d) give the budgets of the four non-zero stresses. The equations 
are nondimensionalized by u T and v. One is struck by the close similarity between the 
pipe and channel results of Mansour, Kim and Moin [22, 23]. Some differences can be 
seen in the budgets of u' r u' r and u' 6 u' e but those were found to be very sensitive to the 
statistical sampling. Many more samples would have to be collected to attribute the 
differences to geometry. In this low Reynolds number range, nondimensionalization 
with respect to u T and v will not completely remove Reynolds number effects which 
probably accounts for most of the discrepancies. 

The budget for the u' z u' z component is dominated by production and dissipation. 
Production reaches a maximum at y + 12 and is balanced by dissipation, turbulent 
transport and viscous diffusion; the minima in turbulent transport and viscous diffu- 
sion occur on each side of the peak in production, or respectively y + ~ 14 and y + ~ 9, 
indicating transport away from the point of peak production. At the wall, production, 
turbulent transport and the velocity pressure-gradient terms are zero by the no-slip 
condition and viscous diffusion balances dissipation. The velocity pressure-gradient 
term is negative throughout the pipe indicating a transfer of energy from u' z u' z to 
the other components; for the streamwise and azimuthal components, the velocity 
pressure-gradient and pressure-strain terms are identical. In the budgets of u' T u' r and 
u' 9 u' e , the velocity pressure-gradient term is the major source of production; with only 
dU z jdr non-zero, production is identically zero for these two components. Figure 3.26 
shows the pressure-strain terms for the normal stresses plotted together. The transfer 
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Figure 3.26: Pressure-strain correlation tensor (II5,j) for the pipe; n5Vr 

ns es — , ns 22 . 


of energy from the streamwise to the cross-stream components is evident. However, 
near the wall ( y + < 12) there is a transfer of energy from the radial component to 
the other two components. This was also observed by Moin and Kim [24] and Moser 
and Moin [25] in channel flows. Moin and Kim attributed this to the impingement of 
fluid coming from the center, hitting the wall, and transferring energy to the other 
components. 

Other terms of the budgets for the u' T u' T and itgu'g components are balanced by 
dissipation, with the remaining terms playing little role, expect in the case of the 
azimuthal component where viscous diffusion balances dissipation at the wall. The 
u' r u' z budget is almost completely dominated by production and the velocity pressure- 
gradient terms balancing each other; for y + > 15, the other terms are virtually zero 
except for a small residual part coming from the turbulent transport term. 



Figure 3.27: Budgets of the turbulent kinetic energy equation (a) and dissipation 
equation (b). B-splines , channel (Mansour et al. [22, 23]) . 
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Figure 3.27 give the budgets for the turbulent kinetic energy (a) and dissipa- 
tion equation (b) ( k and e budgets). Again, similarity with the channel budgets of 
Mansour, Kim and Moin [22, 23] is evident; differences are most probably caused 
by Reynolds number effects. The k budget closely resembles the u' z u' z budget since 
this term dominates u' T u' r and u' 9 u' e . In the k budget however, the velocity pressure 
gradient reduces to the pressure-diffusion term, i.e. II , j = II since the trace of 
the pressure-strain term is identically zero. 

If one defines the dissipation of turbulent kinetic energy, e — e,;/2, a transport 
equation for the dissipation is given by 

jr.e = P] + Pi + Pi + Pi + T e + U e + D c -T (3.13) 


where 


II 
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Mixed production 

Production by mean velocity gradient 

Gradient production 

Turbulent production 

Turbulent transport 

Pressure transport 

Viscous diffusion 

Dissipation 


where 5,-j = 1/2(77, j + Uj,i) is the mean strain rate tensor. 

The e budget also shows great similarity with the channel. Beyond the buffer 
layer, dissipation balances turbulent production; closer to the wall, most terms have 
positive contribution to the e balance with the dissipation term contributing for the 
losses. It should be noted that the dissipation term T is computed by summing all 
the terms in the e equation and not by direct computations (this term involves third 
order tensors which makes its computation non-trivial). Away from the wall ( y + > 30) 
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turbulent production ( P '*) and dissipation (T) dominate the balance. However closer 
to the wall, both P £ and P 2 are of the same order as P*. All terms become zero at 
the wall except for viscous diffusion which balances dissipation. P ’ £ 3 , n £ + T e and D e 
play virtually no role in the balance, except for D £ which reaches a maximum at the 
wall to balance T. 

3.2.6 New structure tensors 

Kassinos, Reynolds and Rogers [16] (hereinafter referred to as KRR) (see also Reynolds 
[35]) introduced a new family of one-point tensors to provide information about the 
structure, or dimensionality of turbulence. Standard turbulence models based on the 
Reynolds stress transport equation (see above) provide information about the com - 
ponentality of turbulence (e.g. R rr = R r $ = R zz = 0 if u' r = 0) but none about its 
dimensionality. 

KRR defined a family of new tensors based on the turbulent stream function; let 

u\ = e i]k4>k,ji with = 0, which implies ip' i kk = -w- (3.14) 

where e ijk is the permutation tensor. Because ipi (dropping the primes for conve- 
nience) satisfies a Poisson equation, it carries non-local information about the velocity 
field (like pressure). Using 3.14, KRR defines 

Rij + Dij + F{j + Cij + Cji — Sijq 2 (3.15) 

where S v is the Kronecker delta, q 2 = Ru and 

R^ = v~v~ = eipqtjtsdwKt Reynolds stress tensor 

Structure dimensionality tensor 
Fij = Structure circulicity tensor 

C^ = d>i,kTpkj Structure inhomogeneity tensor 

All of the above tensors are symmetric, except for Cij. Although 3.15 involves a 
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tensor given by + Cj, which is itself symmetric. 

To compute the stream function ipi, instead of solving directly the Poisson equa- 
tions shown in 3.14, use is made of the stream functions defined in talbe B.l. How- 
ever, if the definitions shown above are to hold, the stream functions have to be made 
solenoidal; let 

, (3.16) 

where 4', are the stream functions from table B.l and 4> , represents the gradient 
of an arbitrary scalar potential. Because e, _,*<!>, = 0 (V x V4> = 0), the velocity 
is identically given by u, = €ijki/>k,j = tijk'&kj- Taking the divergence of 3.16, and 
assuming 0, is constant at the wall, the following Poisson equation and associated 
boundary condition are obtained 


3$ 

$, «. = Wlth ~a~ 

ar 


R 2 


= ~- f 

R> Jo 


1 r R i d’tyr 


dr 


r dr 


(3.17) 


Thus, solving for $ requires the solution of only one Poisson equation instead of three 
for ip. 


Figure 3.28 give the normalized structure circulicity and dimensionality tensors. 
The circulicity tensor provides information about the vorticity field (for homogeneous 
turbulence, KRR show that Fij is directly related to the vorticity spectrum flij ) . As 
its name imply, Dij reveals information about the dimensionality of the flow. As 
shown in figure 3.28 (b), the streamwise ( d zz ) value reaches a minimum at y + fis 10 
indicating that the streamwise motion is nearly 2D, aligned with the 2 axis; circulicity 
(fzz) is also relatively small indicating very little circulation along the axis. This, and 
the peak in the streamwise component of the normalized Reynolds stress (figure 3.29 
(b)) are consistent with the presence of wall streaks; KRR also noted this behavior 
for the channel. They refer to this type of motion as jetal , or jet-like behavior. Also, 
d rr reaches a maximum and d gg a minimum at the same location, indicating that 
the structures vary more rapidly in the radial than in the azimuthal direction, which 
means the streaks have a larger spanwise than radial extent. As seen in figure 3.29 
(a), the wall is also the zone of greatest inhomogeneity, with the normal components 



Figure 3 . 28 : Structure circulicity (/ :J = Fij/F kk ) and structure dimensionality 

( dij = D^/Dkk) tensors; symbols are for radial component , azimuthal , 

streamwise and cross-term 





Figure 3.29: Symmetrized structure inhomogeneity (c i3 + Cji = (C,j + C, ;)/£>*.* ) and 
Reynolds stress (r^ = R i3 /q 2 ) tensors; symbols as for fig. 3.28. 
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of the symmetrized inhomogeneity tensor reaching their maxima at the wall. 

Close to the center, the flow is more 3D with d zz begin somewhat lower than d rr 
and d$$] the structure circulicity tensor also exhibits a similar behavior. The flow is 
still jet-like but less so than near the wall. KRR noted the same behavior of and 
fij for the channel. However an important difference is seen in the inhomogeneity 
tensor: KRR observed that the flow was nearly homogeneous near the center (or at 
least in the log region) which is not the case for the pipe. Both c rr and cge show 
significant inhomogeneity throughout the pipe. 




Chapter 4 
Flow Topology 


This chapter presents a brief description of a topological method used to analyze flow 
structures. This method is then applied to flow fields of the fully developed pipe 
simulation presented in chapter 3. 


4.1 Topological Approach 

One of the fundamental difficulties in analyzing the structure of turbulent flows com- 
puted via numerical simulations is the large amount of information produced by the 
simulations and the apparent lack of coherence of the flow field. DNS typically contain 
several million grid points (2.2 million for the present work), on which three velocity 
components and pressure are computed. To analyze the structure of turbulent flows, 
Perry and Chong [33] introduced a method based on critical point theory. Their ap- 
proach was later extended by Chong, Perry and Cantwell [9] who introduced a general 
classification of flow fields by analyzing the invariants of the velocity gradient tensor. 
This method reduces complex three-dimensional, incompressible flow fields to two- 
dimensional plots of joint probability density functions (PDF) of the invariants of the 
velocity gradient tensor. The information provided by the PDF is essentially qualita- 
tive in nature, revealing general, overall trends rather than in depth, or detailed flow 
features. 
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4.1.1 Description 


A complete description of critical point theory was given in Chong et al. [33] and is 
not repeated here; for the present purpose, a short summary is included for complete- 
ness, and closely follows Sondergaard [42] and Blackburn, Mansour and Cantwell [4] 
(hereinafter referred to as BMC). Define 

A{j = u ij (4.1) 

to be the velocity gradient tensor; again, indices after the comma refer to covariant 
derivative. A{j can be split into a symmetric and skew-symmetric parts 

= Sij + Wij (4-2) 

where 

Sij = \ (Uij + Uj,i) (4.3) 

is the symmetric strain-rate tensor and 

= 2 — u Li) (4-4) 

is the skew-symmetric rotation-rate tensor. The eigenvalues of Aij satisfy the char- 
acteristic equation given by 


A 3 + P\ 2 + QX + R = 0 (4.5) 

where P, Q and R are the invariants of A- VJ and are given by 

P = -Sa (4.6) 

Q = j (j> 2 - 5*5* - (4.7) 

R = ( (-P 3 + 3PQ - S„S, t S it - (4.8) 

* 
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where P, Q and R are respectively the trace, the sum of the determinants of the 
cofactors, and the determinant of A. For incompressible flow, the first invariant is 
identically zero by continuity, 

P = 0 (4.9) 

and only Q and R determine the eigenvalues of A; the discriminant of 4.5 becomes 

D = ?Lb? + qz (4.10) 

4 

D > 0 yields one real and two complex conjugate eigenvalues; D < 0 yields three 
real eigenvalues, and D — 0 gives three real eigenvalues, of which two are repeated. 
For this last case the invariants lie on the lines R = ±(2-\/3/9)(— Q) 3 ^ 2 - 

Figure 4.1 shows the flow topology as a function of Q, R and the discriminant D. 
Each of the four region can be characterized as follows: 

• Stable focus/stretching (R < 0, D > 0): for this topology, the flow spirals in 
toward the origin and flows out along an axis perpendicular to the spiraling 
plane (analogous to vortex stretching). 

• Unstable focus/compressing (R > 0, D > 0): here the flow approaches the origin 
from one direction and spirals out in a plane perpendicular to the approaching 
direction. 

• Stable node/saddle/saddle (R < 0, D < 0): the flow approaches the origin 
from two directions and flows out from the third one (analogous to a stagnation 
point). 

• Unstable node/saddle-saddle (i? > 0, D < 0): the flow approaches the origin 
from one direction and flows out from the other two. 

Other quantities of interest are the invariants of the strain-rate tensor Si } : 

Qs = 

R. = AsijSjkSki 


(4.11) 

(4.12) 
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Q 



e 


Figure 4.1: Invariant space for incompressible flow. The streamlines give the 

flow classification: upper left, stable focus/stretching; upper right, unstable 

focus/compressing; lower left, stable node/saddle/saddle; lower right, unstable 
node/saddle / saddle. 

Because Sij is symmetric, all eigenvalues fall below the line D = 0 in figure 4.1. 
Blackburn et al . showed that joint PDF plots in Q s , R s space provide information 
about the principal direction of strains. Using scaling arguments proposed by Chen 
et al [8], they suggest that fine scale motion should lie far from the origin in Q s , R s 
space. Also, one should note that Q s is proportional to the mechanical dissipation 
of kinetic energy (f> = 2 i/SijSji = —4 vQ s , so that regions of large dissipation are 
identified by regions of large negative values of Q s . 

Also, plots of Q s versus the second invariant of W{j , Q w , 

<J» = -jl VijWji 


( 4 . 13 ) 
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reveal relative importance of straining versus rotation. Since Q w is proportional to 
the enstrophy (ujju;,*), it provides information about the rotation of the flow. 


4.1.2 Results 

The method of using invariants of the velocity gradient and other associated tensors 
was used by Soria et al. [43] to study the topology of incompressible mixing layers, 
Chen et al. [8] for compressible and incompressible mixing layers, Sondergaard [42] 
for transitional wakes, and BMC for fully developed turbulent channel flow. 

Figure 4.2 (a)-(d) give the joint PDF in (Q,R) space for a single realization at 
four radial locations, each in a different region of the flow: viscous, buffer, logarithmic 
and wake region. For all regions, except close to the wall, the distribution in ( Q,R ) 
space shows a preference for the second and fourth quadrant, while close to the wall 
the distribution is almost circular with no preference for any quadrant. Cantwell 
[6] showed that under certain assumptions (one of them being inviscid flow), the 
preference of the invariants to line-up in the second and fourth quadrants in a “skewed 
teardrop” shape is an actual solution to the Navier-Stokes equations. These results are 
in agreement with the findings of BMC in channel flow. The high level of noise near 
the center is again due to the fewer number of modes which degrades sampling (c.f. 
figure 3.10). The behavior near the center is consistent with the topology observed 
in free shear flows, i.e. stable focus/stretching and unstable node/saddle/saddle. 

Figure 4.3 (a)-(d) show plots of the second and third invariants of the strain-rate 
tensor, for the same four radial locations. Since the strain-rate tensor is symmetric, all 
the eigenvalues are real and lie in the region of negative discriminant (D < 0). Away 
from the wall there is a clear preference for the unstable node/saddle/saddle topol- 
ogy. As the wall is approached, the preference for the unstable node/saddle/saddle 
topology is no longer observed, and very close to the wall the probability tends to 
cluster close to the R s = 0 line. This is in agreement with BMC’s findings for channel 
flow; they also showed that at the wall R s = 0. The R s = 0 line corresponds to the 
case of purely two-dimensional strain; if one orders the eigenvalues of Sij such that 
Aj > A 2 > A 3 , then for R s — 0, correspond A 2 = 0. As with the structure tensors 
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(a) r= 169.1 r/D = 0.055 



(b) y*= 122 r/D= 0.1789 



Figure 4.2: See caption page 85. 






4.1. TOPOLOGICAL APPROACH 


85 



R 



R 


Figure 4.2: Contour plot of joint PDF (log 10 ) of Q and R. Contour levels are from 
the exterior toward the interior: (a) and (b) 0.5-2. 5, (c) 0. 5-3.0 and (d) 0.5-3. 5 by 
half decade increments (0.5). 
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(a) r = 169.1 r/D = 0.055 



(b) y* = 122 r/D= 0.1789 



Figure 4.3: See caption page 87. 








Figure 4.3: Contour plot of joint PDF (log 10 ) of Q s and R s . Contour levels are from 
the bottom toward the origin: (a)-(c) 0.5-2.5, (d) 0.5-2.0 by half decade increments 
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(section 3.2.6), the invariants provide information about the dimensionality of the 
flow field. 

Plots of SijSij versus —WijW{ 3 are shown in figures 4.4 (a)-(d). Very close to 
the wall, the results collapse onto a single 45° line, indicating that dissipation and 
enstrophy are identical. This type of topology is analogous to a vortex sheet. In 
the buffer layer, the results still maintain a 45°-like behavior but much more losely, 
whereas further out in the flow, the results look like that of free shear flows (e.g. 
Sondergaard [42]). Plots (a) and (b) also show an interesting feature: the contours 
never reach a value of zero strain, but can reach zero rotation. This implies that 
although pure strain is possible, no point in the flow is in solid-body rotation. 


4.2 Vortices 

One of the more interesting features in pipe flow and wall-bounded flows in general, 
are streamwise vortices. There lies a great difficulty in identifying vortices in turbulent 
flows, since there is no universally adopted definition of a vortex. For example, if one 
relies solely on vorticity magnitude, the region near the wall possesses the largest 
vorticity but yet has no vortex, hence the difference between vortex-sheet and tube- 
like vortices or vortical eddies (see for example Chong et al. or Jeong and Hussain 

[14])- 

A topological method, such as the one used in section 4.1, can be used to identify 
the regions of the flow where tube-like vortices are present. Figure 4.1 shows that 
regions with positive discriminant have streamlines spiraling around the principal 
axis. So regions of the flow with D > 0, i.e. where the eigenvalues of Aij are complex, 
should correspond to regions of (tubular) vorticity. However, Jeong and Hussain [14] 
suggested that this criterion (D > 0) was too general indicating that streamlines can 
spiral without being a vortex. They proposed an alternate approach based on the 
eigenvalues of the symmetric tensor 


Tij = S ik S ki + WikWkj 


(4.14) 
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which represents the nonlinear source term of the transport equation for S{j. They 
define a vortex as a region with two negative eigenvalues of ordering the eigen- 
values such that Ai > A 2 > A 3 , the definition is equivalent to requiring A 2 < 0. The 
potential problem with Jeong and Hussain’s definition lies in the arbitrariness of the 
definition of Tp. A comparison of both approaches (1) D > 0 and (2) A 2 < 0 is in 
order. 

Figures 4.5 (a) to (c) respectively show iso-surface plots of D, A 2 and streamwise 
vorticity u z . The direction of the flow is from left to right. Contours levels were 
adjusted to get the best possible agreement with the vorticity plot. Clearly, the two 
approaches give similar results when comparing with the vorticity. Most of the dom- 
inant features are captured by both methods. As Jeong and Hussain point out, when 
compared with the plot of D , the eigenvalue plot does leave out some information 
since it was designed to be more restrictive than the discriminant. Both topological 
approaches do not include the region of high vorticity at the wall (vortex sheet), such 
that all features revealed by figures 4.5 (a) and (b) are away from the wall, as de- 
sired. The vorticity plot show tube-like structure oriented with the flow, which tend 
to extend in the flow away from the wall. This behavior is also captured by both 
topological approaches. 

So, both approaches seem to paint a good overall picture of the organized struc- 
tures of the flow, although they both require some adjustment of the contour levels 
in order to get a good match. Within the adopted contour levels, the discriminant 
method gives marginally better results, matching the actual vorticity contours more 
closely, than Jeong and Hussain’s method. 

Wall vortex 

In section 3.2.2 (see figure 3.14), a short, high intensity wall vortex was shown to 
produce a high speed streak as it moved downstream in the pipe. A plot of the 
discriminant for this feature in given in figure 4.7 (a). The discriminant reaches a 
maximum at the location where the vortex is known to reside (see above). Chacin 
and Cantwell (private communication) observed by studying the turbulent structure 
in boundary layer, that the discriminant tends to peak where turbulent production 
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(b) A 2 



(c) w 2 


Figure 4.5: Iso-surface plots of: (a) the discriminant Z), (b) the intermediate eigen- 
value A 2 , and (c) the streamwise vorticity u z . 
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Figure 4.6: Maximum value of the discriminant versus y + for the region containing 
the wall vortex of figure 3.14. 


(Pk) also reaches a maximum. Figure 4.7 (b) shows a discriminant plot for y + ~ 12.5, 
where turbulent production is maximum (see figure 3.27), and the overall trend of 
the maximum value of the discriminant is shown on figure 4.6, clearly showing the 
peak at y + ~ 12.5. The maximum value is taken from the region containing the wall 
vortex. The results are consistent with Chacin and Cantwell’s observations, with a 
maximum discriminant an order of magnitude larger than that observed at y + 5. 
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Chapter 5 
Conclusions 


This chapter concludes the present work with an overview of the development of the 
numerical method and the results of the pipe flow simulation. Some recommendations 
for future work are also made. 


5.1 Numerical Method 

A new numerical method for the computation of incompressible flows in cylindri- 
cal geometries was developed. The method is based on a vector Galerkin weighted 
residual method, using divergence-free expansion and weight vectors. With conti- 
nuity satisfied a priori and pressure dropping out of the formulation, this approach 
reduces the number of unknowns from four to two, making memory management 
straightforward. The method makes use of b-splines polynomials to represent the 
radial direction and of Fourier transforms for the other two directions. This approach 
was found to possess spectral-like accuracy while providing a flexibility unavailable 
in standard spectral methods. In fact, comparing with results obtained from a finite 
difference simulation demonstrated the superiority of the present approach. The flex- 
ibility of the method was used to guarantee that the expansions were regular near 
the origin by satisfying a set of conditions. Also, a procedure was implemented to 
alleviate the restriction put on the time step by removing unnecessary b-splines near 
the origin. 
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5.2 Pipe Flow Simulation 

A direct numerical simulation of incompressible pipe flow was successfully carried out 
using the method described above. The simulation was run at a Reynolds number of 
5600 based on the diameter and bulk velocity, in line with several experiments and 
other simulations of both pipe and channel flows. The most important observation 
is that pipe flow and channel flow have very similar turbulent statistics. Differences 
between both flows are nevertheless apparent in the mean profiles with pipe flow not 
obeying the log-law (in this low Reynolds number range). The spectra and two-point 
correlations revealed that the flow was well resolved, although the pipe length is a 
little too short. However this does not seem to have played significantly as all results 
agree very well with experiments. 

Turbulent intensities and higher order moments agreed generally well with ex- 
periments and other simulations. For higher order moments however, the present 
simulation with its high resolution and accuracy revealed some shortcomings of near 
wall measurements. Furthermore, it was shown that flatness levels predicted with the 
b-splines approach are more reliable than the ones experimentally measured or com- 
puted by finite differences; high levels of flatness in the near wall region are physical. 
A possible explanation for the formation of high speed streaks in the near wall region 
was provided by the presence of short streamwise eddies (/* 65 wall units) inducing 

positive streamwise velocity perturbations in its wake. Budgets of turbulence trans- 
port equations were again very similar to the channel’s, which means that standard 
turbulence models should perform reasonably well in pipe flow, having proven their 
adequacy in channel flow. New structure tensors were computed and revealed that 
contrary to channel flow, the logarithmic region of pipe flow is not homogeneous (for 
this Reynolds number). 

A topological method consisting in a classification of the second and third invari- 
ants of the velocity gradient tensor revealed that the flow (away from the wall) tends 
to adopt a preference for the stable focus/stretching and unstable node/saddle/saddle 
topologies, similar to that of free shear flows. Closer to the wall, the flow was found 
to have no preference. It was also shown that iso-surface plots of determinant proved 
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very useful in visualizing flow structures, mainly vortical eddies. 


5.3 Recommendations for Future Work 

Building on the present work, several avenues could be taken in the future: 

• Study the time evolution of turbulent structures in pipe flow. Such work 
would clarify the relationship (causality) between streamwise vorticity and wall 
streaks, while enlarging the statistical sampling of the simulation, thus resolving 
some of the observed irregularities of the present simulation. 

• Increase the Reynolds number of the pipe. This would provide for an inter- 
esting comparison between two different Reynolds numbers, and should be less 
sensitive to the present pipe length. 

• Explore several classes of free shear flows, mainly jet flow. By implementing a 
set of free shear boundary conditions and using the pipe flow simulation as a 
starting point, it will be possible to study the development of a fully turbulent 
round jet. Free shear flows should also not be affected by the short domain size, 
since streamwise coherence should not by preserved in the absence of the wall. 

Some work could also be performed on the computer code itself, such as improving 
the performance of the parallel version. 




Appendix A 

Few Facts About B-splines 


This appendix introduces some of the basic properties of b-splines. The information 
provided here should be sufficient to understand most of the issues relevant to this 
work. For more information on the subject, the reader is invited to consult de Boor 
[10] (primarily) and also Shariff and Moser [41]. 

A.l Background 

Although similar in concept to standard finite element polynomials, b-splines offer 
spectral-like accuracy and are C 5-1 continuous, where S is the order (or degree) 
of spline being used; this means that derivatives of velocity, such as vorticity, are 
smoothly and accurately represented. The higher the order, the greater the accuracy, 
storage requirements and computational cost. A trade-off must be reached between 
reasonable accuracy and fast computations. The spectral accuracy of b-splines is 
shown on figure A.l, where both eigenvalues of first and second order derivative op- 
erator are shown. For higher degree splines (third or fourth) over two third of the 
spectrum of the first order derivative operator (or modified wave number), and vir- 
tually the entire spectrum of the second order operator are reproduced accurately by 
the b-splines. Standard second order finite difference results are shown for reference. 

By construction, a b-spline is required to have zero value and zero first 5 — 1 
derivatives at the edges of the intervals where it has support (see figure A. 2); as such, 
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£1 



(b) 


Figure A.l: Spectral accuracy of b-splines of degrees 1 through 4: (a) Modified wave 
number, or eigenvalue of the first order derivative operator; (b) eigenvalue of the 
second order derivative operator. The straight line is the exact or spectral limit and 
higher degree are closer to exact. The dashed line is the second order finite difference. 


b-splines are the smoothest type of splines. As their name imply, b-splines form a 
basis, such that any smooth function / can be represented as a linear combination of 
b-splines; defining a; to be the expansion coefficients and g t the b-splines polynomials, 
we write the linear combination 

/( r ) = Y, a i9,( r ) (A.l) 

/=i 

where N T is the number of b-splines and it is implied that the b-spline basis is built 
such that 

= 1 (A.2) 

(=i 

The at coefficients can then be determined by using an L 2 projection 


a = A *b 


(A. 3) 
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Figure A. 2: Quadratic b-splines on an equispaced grid. The x symbols represent the 
knots which delimit the intervals over which b-splines have support. 


where A is the N r x N r mass matrix 


f L 

A = {a n } = / S ( ,( r )s,(r) r dr 

J o 

(A.4) 

fL 

b = {&('} = J g l ,i r )f ( r ) rdr 

(A.5) 


and a = . ,ajv r }- 


A. 2 B-splines Construction 


Traditionally, splines are constructed by solving a linear set of equations. However, 
in the methodology established by de Boor [10], constructing b-splines is done by 
solving a simple recurrence relationship. First, defining 


1 ,7)1 <r < 7)1 + ! 

0 , otherwise 


(A.6) 


where gi- n is the / th b-spline of degree n and Tji is the knot coordinate. Knots are used 
to delimit the intervals over which b-splines have support. It should be noted that 
A.6 is a direct consequence of A. 2. The recurrence relation becomes (see de Boor 
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page 131) 


9l;n+l(r) - - — gi- n (r) + 


Vl+n- 1 — r 


- ii^ 9, * Ur) (A.7) 

To construct b-splines near the boundaries, de Boor introduced the concept of 

bl“Z ??? knots ' Al the boundaries fcwer continuity const ™" ts ■“» ^ 

t ’ n u u ^ USed iB lhe CO " tei “ ° f the reCurren “ Nation A .6 and A 7 

ilrat t Li d f ° f T 1 "" 17 ° f b ' SP,meS - FigUrC A ' 3 Sh ° WS «“ ™l«Ple 

OT * b ' SP,ine ° f * ‘ here “ 5 -outtipk knots at the 


n 1^2,113 


-X — 

n 4 


-x — 

n 5 


-X — 


% 

r 


* * x 

118 ^ Hio Tli,.Tll2.ni3 


tefj e af,he bo n „°nda U dl in C ° nS ‘ rUCtin * b -is of figure A.2. Note the degenerate 


in 


fit" the Pr ° CedUre ' let " COnStr “ Ct «ret two b-splines shown 
* " ° rder l ° PreVent divis ™ s ^ zero ’ the first step consists in letting 

S , :0=feO = ° (A. 8) 


9l;l = 0 


(A. 9) 


The zeroth order or constant b-splines are shown in figure A 4 (al with ,h 

“■* ■" - * - — «. w «. l 

rom the zeroth order one using A.6 and A.7. They are given by 


02 , -i (r) = 



93;l(r) = J 


Vi-r)3 

V5~r 

V5-V4 

0 


, V3 < r < r) 4 

, otherwise 

(A. 10) 

> 93 < r < r )4 
i V4 < r < t) 5 
, otherwise 

(A. 11) 
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Figure A. 4 (b) shows the first linear b-splines. It is interesting to note that lin- 
ear b-splines are identical to standard linear finite element polynomials. Lastly, the 
quadratic b-splines are obtained from the linear b-splines; the functions are given by: 


9i*(r) 


{v 4-r) 2 
(r) 4“t?3) 2 

0 


, 7]3 < r < T )4 

, otherwise 


(A. 12) 


9 2 : 2 (r) 


' (r-V3)(Vi-r) , ( t ? 5 -r)(r-rj 3 ) 

(V4-V3) 2 (V5-r)3)(V4-V3) 

< (vs-r ) 2 

{Bs-V3){Vb-V4) 


{ 0 


,T) 3 <r<T) 4 
,Tj 4 <r <r) 5 
, otherwise 


(A.13) 


with the functions shown in figure A. 4 (c). 


A. 3 Support Rules 

Because b-splines have local support, the bandwidth of the mass matrix, or any other 
matrix which results from an inner product operation, will depend on the degree of 
b-splines used. Figure A. 2 shows quadratic b-splines: each b-spline has support on 
three intervals, which are delimited by knots. For example, b-spline 6 has support on 
intervals shared by splines 4, 5, 7 and 8. This implies that the inner product computed 
in A. 4 will have a bandwidth of 5. In general, a b-spline of degree S will have support 
on S + 1 intervals (except near the boundaries where b-splines have support on fewer 
intervals since there are fewer continuity constraints to be met) translating in the 
inner product having a bandwidth of 25 + 1; the total number of non-degenerate (or 
single) knots, JV*, is then Nk = N r — 5 + 1, where N r is the number of b-splines. 

The behavior of the b-splines near the boundaries of the domain is extremely im- 
portant in understanding how to manipulate the expansion coefficients to implement 
the boundary conditions. A general formula for the b-spline support at r = 0 is given 
by, 

5 ^( 0 ) = 0 for / > q + 1 and q = 0, 1, 2, . . . , S (A. 14) 

where the superscript q refers to the <? th derivative. Similarly, at the other boundary, 
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r = 1 

( 7 , (,) ( 1 ) = 0 for l>N r + q and q = 0, 1, 2, . . . , S (A.15) 

Both A. 14 and A.15 can be tabulated to yield: 

Table A.l: Support rule for b-splines at r = 0. For r = 1 the table is identical except 
that the index i in is replaced by N r - (i- 1). The “x” represents non-zero values. 


q 

si 9) 



*?> 


>) 

9s 

9s + 1 

to 


0 

X 

0 








1 

X 

X 

0 







2 

X 

X 

X 

0 






; 




* . 






5 

X 

X 

X 

X 


X 

X 

0 



Using A. 2, A. 14 and A.15 two important special cases can be derived; first 


Si(0) = 9N r ( 1) = 1 


(A- 16) 


and taking the derivative of A. 2 to yield g\{r ) = 0, we get 


fifi(O) = -£ 2 ( 0 ) 

5N r (l) = -5AT r -l(l) 


(A- 17) 
(A- 18) 





Appendix B 


Regularity Conditions 


In section 2.3.2 of chapter 2 we introduced the concept that any vector field written in 
cylindrical coordinates must meet a certain number of conditions in order to be regular 
near the origin. The purpose of this appendix is to examine in greater details how 
these conditions come into play when the vector expansion functions are constructed. 
An example is also included to delineate the whole procedure. 


B.l Regularity Conditions 

Shariff (private communication) derived the regularity conditions presented in this 
appendix. For completeness, we include the derivation of the conditions for u z ; the 
conditions for the other two components can be obtained in a similar manner. Let u 
be any vector field in cylindrical coordinates: 



(B.l) 

u T = u r (r-,kg,k z )e l( ' ke6+kzZ ^ 

(B.2) 


(B.3) 


Shariff’s method consists in starting from a vector field in cylindrical coordinates and 
transforming it to cartesian coordinates by paying attention to the behavior at the 
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origin. Define 


A T> 

u z ~ ar 


(B.4) 


where a is a complex constant which depends on p, kg and k z . In the complex 
plane, the transformation between cartesian and cylindrical coordinates is given by 
x = r cos 6 and iy — ir sin 0 such that e' keS can be rewritten as 


e »M _ £ 
r 

where ( = x + iy and r = |£| = y/x ^ + y 1 . u z becomes 


(B.5) 


it, = 


(B.6) 


Clearly, (^ k ° is regular for kg > 0 and r p ~ ke is regular only for p — kg = 0, 2, 4, 6, . . . 
since derivatives of odd powers of r can involve terms in r -1 / 2 which are not regular 
at r = 0. So, we get 


p=kg+2n where n = 0,1,2,... (B-7) 

and the regularity condition for u z follows 

ii z (r;kg,k z ) = a(k 9 , k z ) r ke P z (r 2 -,kg,k z ), kg > 0 (B.8) 

where P z (r 2 ; kg , k z ) is a polynomial in r 2 . Regularity conditions for the other compo- 
nents are given by 


u r (r-kg,k z ) = b(kg,k z ) r ke 1 P r (r 2 - kg,k z ), kg > 1 (B.9) 

u e (r;kg,k z ) = ib(kg, k z ) r ke ~ l Pg(r 2 -, kg, k z ), kg > 1 (B.10) 

and when kg = 0, 

u r (r; 0, k z ) = 6(0, k z ) rP r (r 2 : 0, k z ) (B. 1 1) 

ug(r-,0,k z ) = c(k z ) rPg(r 2 ;0,k z ) (B.12) 
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where a, b and c are constants that can depend on the wave numbers, Pi(r 2 ;k e ,k z ) 
are polynomials in r 2 , and P,(0; k z ) = 1. Conditions B.8 through B.12 do not 
account for negative azimuthal wave numbers, since these are obtained by symmetry 
and not by direct computations. 


B.2 Regularity of the velocity vectors 

After dividing the expansions into plus and minus modes, the following vectors were 
adopted 


V x 

'Srf(r-,k e ,k z ) 

(B.13) 


l 0 \ 


V x 

0 

(B.14) 


—k z rg t 

and 

u i~(r;k$,k z ) - V x ^f(r;ke,k z ) 

( -iff, ) 

9i 

0 / 

and because the above vectors are incomplete when k z = 0 * 

u; + (r; k$,0) = Vx$f(r;l ( ,0) 

/ 0 \ 


V x 


0 

V - r 9> ) 


(B.15) 

(B.16) 


(B.17) 

(B-18) 


‘Since k z = 0 causes completeness problems, one might question the need for k z to be present in 
B.14 in the first place. With k z present in B.14, and later in B.21, the mass and viscous matrices are 
independent of the sign of k z (they depend on some powers of k z 2 ), which reduces the computations 
since both matrices only have to be computed for half of the k z spectrum. It should be noted that 
the exact solution to the Stokes problem is also independent of the sign of k z , a fact which must be 
reflected in the numerical approximation. 
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and u/ unchanged. 

In order for the velocities u/ + and U/“ to be regular, the vector stream functions 
i&l and must also be regular, since taking the curl will not affect regularity 
of the vectors either favorably or unfavorably. It is clear that the first condition 
which requires the radial and azimuthal components to be constrained to each other 
is satisfied (B.9 and B.10). The second condition which requires a certain behavior 
in r near the origin is generally not satisfied by the b-splines (g t ). To alleviate this, 
the expansion coefficients af ml are constrained such that the linear combination of 
b-splines does have the correct behavior in r. For the vectors B.14, B.16 and B.18 
this constraint can be written as 

Yl a fmi9,(r) ~ r k °- l P(r 2 ) (B.19) 

/=! 

To ease implementation of the constraint in the computer code, special care was 
taken when designing the expansion vectors so the constraint would be identical for 
both the plus and minus modes. B.19 implies that not only the values of the splines 
must be constrained, but also a certain number of derivatives. If S is the order of 
b-splines being used, knowing that there are 5 — 1 continuous derivatives at any given 
point, B.19 can be generalized as follows 

Z,af m M0) = af„M0) = cf mI ~ r k '~'P( r J )| _ o 

(0) + af m2 si(0) ~ P(r J ))| r=o 

= <i*r"(0) + • ■ • + af mS g ( t'\0) ~ J£-(r t »-'/>(r*))| r _ o 

r (B.20) 

where A. 14 was used to simplify the left hand side 

When the right hand side of B.20 is nonzero the af m i coefficients are unconstrained, 
otherwise the sum is set to zero. From B.20, when k$ > 5+1, the first 5 coefficients 
are automatically zero. The constraints are implemented by modifying the necessary 
lines of the mass and viscous matrices in 2.18 (at most the first 5 lines), such that 
the modified mass and viscous matrices reflect the constraint in B.20. 
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For the cases when kg = 0, the constraints in B.20 cannot be applied to B.14, 
B.16 and B.18 since they involve terms in r _1 . Two additional vectors must be used 
for this case 

( 0 \ 


u/ + = 


k z 9i 

\ 0 } 


(B-21) 


with u;~ still unchanged. When k z = 0 the vectors are again incomplete and another 
set of vectors is used 


/o\ 

1 i 

0 ^ 

9i 

II 

1 

P 

a 

0 

\oJ 

f i 

\?) 


(B.22) 


For these last cases (kg = 0), the expansion coefficients aQ mi are constrained 
by using rP(r 2 ) as the right hand side in B.20 instead of r* e-1 P(r 2 ). Again, the 
constraints are identical for both the plus and minus modes. 

With the expansion vectors given in table B.l and the constraints B.20, we observe 
that constraining the first 5—1 derivatives also results in constraining the first 5 — 1 
derivatives of u r and ug, but only the first 5 — 2 derivatives of u z , since u z already 
involves a derivative of g r Even though ug also involves a derivative of g n it is 
multiplied by r, such that at the origin this term vanishes. 


Mode Redundancy 

A problem exists with the expansion vectors shown in table B.l for kg > 0: there 
are four coefficients which have support at the boundary, but only three velocity 
components. This was found to cause conditioning problems with the mass matrix 
rendering it singular (uninvertable). To alleviate this problem, we simply set 

a JmNr =0 f° r kg > 0 (B.23) 

A possible explanation for this problem comes from the presence of derivatives of 
b-splines in both vectors. Recall the expression for the number of knots in section A 


N k = Nr - S + 1 


(B.24) 
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Table B.l: Summary of expansion vectors and associated stream functions 



Since the number of knots is constant, with derivatives of b-splines of degree 5 exactly 
representable by b-splines of degree 5 — 1, in order to keep the number of knots 
constant, only ./V r — 1 b-splines of degree 5—1 are necessary, hence the extra b-spline. 


B.3 Regularity of the weight vectors 

Because a Galerkin method is used, the weight vectors w/< are the complex conjugate 
of the velocity vectors u; given in table B.l; this means that regularity of the weight 
vectors must also be ensured. Here again, the idea is to take a linear combination of 
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weight vectors, similar to B.19, such that the linear combination be regular 

Pnl'^n (B.25) 

n— 1 

where is the regular weight vector. With the weight vectors obtained from the 
velocity vectors, B.25 is equivalent to B.19 and B.20 such that 

E/Wr)-^- 1 /^ 2 ) (B.26) 

71=1 

except when kg = 0 where the right hand side of B.26 becomes rP(r 2 ). By expanding 
B.26 as in B.20, a system of the form G(3 = 0 is obtained, where (3 is a matrix 
whose columns are the null vectors of G, the matrix of constraints. Each component 
of the null vectors are the /3 coefficients in B.26. Because the null space of G is 
non unique, the following choices are made depending on the azimuthal wave number 
and the degree of b-splines being used (those /3 coefficients are computed once at the 
beginning of the code). 

•kg and S odd: 

1 0 0 

A i o o 

g\ g' 2 0 0 0 ... 0 0\ 0 1 ... 0 

s'" < 93" <' 0 ••• 0 0 A* fti 0 

; : 000 

9? 9 g? st gf ... & 0) • 

Ad A(d-i) • • • 0 
^00 1 

where p = S — 2, d = (S — l)/2 and all the splines are evaluated at the origin, i.e. 
g\ n ^ = | . With this particular form of null space, B.25 is written as 

wfi = w/A + E P(iLti){ i±xziL ) w t l' < S and /' odd (B.28) 

n=l'+i 2 2 

n even 
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•kg odd and S even: 


/ 1 


9\ 

92 

0 

0 

0 . 

.. 0 \ 

g\ r 

S? 


97 

0 . 

.. 0 


9 ( 2 P) 

<4 p) 

9 ( : ] 

<4 P) • 

•• 9 { I ] J 


Pii 

0 

@12 

0 


0 

0 

1 

/?21 

0 


0 \ 
0 
0 
0 


\0 Id 02(d-l) • • • ftd\) 

where p = S — 1 and d = S/2. In this case, is written as in B.28. 


•kg even and S odd: 


/ 0 0 

1 0 


f 9i 

0 

0 

0 

0 . 

.. 0 \ 


Pu 

0 

9’i 

// 

92 

9 3 

0 

0 

.. 0 


0 

1 

\9 i P) 

d p> 


9 { I ] 

5i P) • 

•• 9 ( s P) ) 


fin 

p2l 


0 0 
V ft Id 02(d-l) 


where p = S — 1 and d = (S — l)/2 and B.25 is given by 


5 

W* = w,.* + 52 P(iL)(l±n=lL)™n 

n=:l '+ 1 2 2 

n odd 


0 \ 

0 

0 

... 0 
0 

1 

■ - - Pdlj 


0 (B.29) 


(B.30) 


l' < S and l' even 


(B.31) 
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•kg and S even: 








° 

0 

0\ 







1 

0 

0 

9i 

0 

0 

0 

0 

0 0\ 


Pn 

0 

0 

aH 

92 

9 3 

0 

0 

0 0 


0 

1 

... 0 

* 





* 


012 

021 

0 

wi p) 

9 { 2 P) 

^ P) 

9? 

9? 

••• 9(s- 1) 0^ 











01d 

02(d-l) 

... 0 







0 

0 

1/ 

:e p = 

5-2 

and d 

= (5 

— 2) / 2 and the weigth vectors 



= 0 


(B.32) 


S-2 


W,, 


= W,^ + /' < 5 - 2 and /' even 


n=/'+l 
n odd 


(B.33) 


So, for a given b-spline of degree S , there are two families of 0 coefficients: one 
for even and one for odd azimuthal wave numbers. 


B.4 Example 

The best way to illustrate how to implement the above conditions is through an 
example. Suppose we wish to solve for 


Aa = b (B.34) 

Typically, this problem can arise when solving the Poisson equation for the pressure. 
Consider 5 = 3 (cubic b-splines) and kg = 1; B.27 reduces to 


- 0 


( 9 [ 0 ) 


1 °\ 
Pn 0 
0 l) 


(B.35) 



116 


APPENDIX B. REGULARITY CONDITIONS 


Solving yields j3 u = —g[/g 2 - Consider now the upper portion of the original system 
B.34; writing the matrix A with a bandwidth 25+1, 


f a \\ a\2 G13 014 ^ 


( ai ) 


( bl \ 

a 2 i a 22 a 2 3 a 24 a 2 $ 


«2 


b 2 

O31 032 <*33 O34 035 036 


«3 


b 3 



i ) 


: / 


(B.36) 


The first step is to make w regular by taking a linear combination; using B.28 


W, = Wi + /?11W 2 (B.37) 

W3 = W3 


which means that only the first line of B.34 will be modified, since the index of the 
weight functions (/') corresponds to the line number of the matrices or the vectors. 


f dn 

^12 

&13 

«14 

<215 



/aj\ 


( bl \ 

021 

022 

<223 

«24 

<225 



«2 


b 2 

a 31 

<232 

<233 

a 34 

<235 G 36 



«3 


b 3 




4 . 


/ 


i ) 


: / 


(B.38) 


where following B.37, dij = a X j + /3 n a 2j and 61 = 61 + /?n 6 2 . 

The second step is to ensure regularity of the expansion functions by constraining 
the ai coefficients. Using B.20 the only relevant constraint, i.e. for which the right 
hand side of B.20 is zero, is 

ot\g[ + a 2 g’ 2 = 0 (B.39) 


and building that constraint into the system of equations yields 


( dn 

^12 

013 

di 4 

< 2 l 5 

s[ 

9 2 

0 

0 

0 

031 

032 

033 

034 

<235 <236 




\ 

(<*i\ 


/M 


a 2 


0 




b 3 

/ 

: J 


: / 


(B.40) 


which is the system that will ensure that the expansion possesses the correct behavior 
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at the origin. 

We conclude this section with two important comments. First, from this example 
it is clear that with this particular form of null space, no information was lost by 
imposing the regularity conditions, i.e. the expansion is still complete. Even though 
B.39 was imposed by removing the second line of the matrix, the information that 
was contained was absorbed before hand in the first line of the matrix with B.37. 

Second, the conditioning of the matrix can be adversely affected by the regularity 
conditions. This is especially true if splines of high degree are used, since as can 
be seen from B.20, with higher degree splines it is possible to impose constraints of 
higher order. With derivatives being approximated by 

4 9) ( 0) - 0(Ar _<? ) (B.41) 

it is easy to see that when A r is small, higher order derivatives can quickly become 
very large ^ . The consequence of this is to limit the highest degree of b-splines to 
about four or five. In order to preserve accuracy of the solution, double precision 
must be used throughout the computations; the standard approach on the CRAY 
super computer which consists in truncating variables to single precision when stored 
on disk must be avoided (see Moser and Moin [25]). 


tThis problem is made more difficult by the fact that even when the regularity conditions are not 
applied, the conditioning of the matrices is degraded when A r is made smaller. This is a standard 
result in finite element analysis which is also observed here (see Strang and Fix [46], section 5.2). 




Appendix C 
Implementation 


In this appendix, we present the full definition of the mass and viscous matrices 
and complete the details on how the nonlinear term is computed. Extensive use of 
Mathematical a symbolic manipulation package, greatly reduced the effort in deriving 
what follows. 


C.l Mass and Viscous Matrices 


The mass and viscous matrices are constructed from the following elemental matrices 
which are computed (to machine accuracy) once at the beginning of the code using 
Gauss quadratures. The complete mass and viscous matrices have to be reassembled 
at every time step, since it would require too much memory to store them. The 
following expressions are obtained by substituting the vector expansion functions 
(see table B.l) in 2.19 and 2.21; noting that g, = g t (r), the following are defined: 


mi 


m 3 


m 5 


Jo 

rR2 

= J 9,9 t . r & 

rR2 

=J o 9,%r 2 dr 


m 2 


m 4 


m 6 


r dr 


rtu 

= i 9 ' 91 ' 

r R2 

= / 9 l l 9 l , l r z dr + Rlg l {R 2 )g l ,{R2) 

J 0 


rR 2 

=L 9 ^' ] 


1 dr 
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m 7 = 


m 9 = 


nin — 


mi3= 


fR2„ „ 

l ^rfr-s,'(0)s,,'( 0) 

/•H2 

I 9 " 9 / r dr- R 2 g l ,, (R 2 )g l /(R 2 ) 

J 0 
rR2 

/ s/s,, 

*/ 0 

rR2 


rfl2 / / 

' P/ 5»/ / 

m 8 = / ■■ ar 


=/„ 


r /?2 

mio=jf g"g,,"r 3 dr- 

(Dlg°"(H 2 )g/(R 2 ) + Rh,"(Ri)g„(R 2 )) 

j-R 2 

mi2= / S,S,/ rfr 

JO 

rR2 

m 14 = y o g/g/'r 2 dr - /?^/(/? 2 )^/( J R 2 ) 


furthermore, the following boundary terms are also defined (matrices which only 
have non-zero elements on the last line, i.e. /' = jV r , by virtue of A. 14) 

bt, = 9l (R 2 ) gil (R 2 ), bt 2 = R 2 g/(R 2 )g l ,(R 2 ), bt 3 = 

)s,,(# 2) . /// p /o \ s/(#2)s,,(#2) g,(R2)g,,(R2 

Note that the boundary term in m 7 was obtained by applying L’Hopital rule to 
boundary terms containing a l/r n factor. Other boundary terms at r = 0 cancel 
either because they are absorbed by the imposition of the regularity conditions, or 
contain an r n factor. 



C.l.l Mass Matrices 

The mass matrix is assembled from the following matrices (see 2.24 and 2.25). 

• k z ^ 0 , kg 7 ^ 0 : 

A+ = = k 2 \(kg 2 - l)m 3 + m 4 ] 

AI = {aZ ri } = \2k z 2 m 3 + (kg - l) 2 mi + m 2 + (1 - fc/)bt,] 


At = {at,,J = k 2 (k 8 + l)m 3 + m 5 
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A + = {a + j, ( } = k z 2 + l)ni3 + nig 

• k z = 0, ke 7^ 0: 

A+ = [(A ;* 2 - l)m 3 + m., 

A I = [(^ - l) 2 mi + m 2 + (1 - &o)bti 

At = a; = o 

• k z ^ 0, k$ = 0: 

A+ = Al = A; = k z 2 m3 
A I = [2fc 2 2 m 3 + m 1 + m 2 + bti] 

• A;* = 0, = 0: 

A+ = m 3 
AI = mi 

Al = a; = o 


In order to minimize storage, the mass matrix is stored in banded form, following 
the structure shown on figure C.l. Each element of the matrix takes its value from 
the four sub-matrices defined above. The elements denoted by an “X” appear after 
the imposition of the regularity conditions. As seen in the example of section B.4, 
imposing those conditions will locally increase the bandwidth of the matrix. Because 
of this unusual feature, special banded matrix solver and matrix-vector multiply 
routines were written to account for these extra terms. With the mass and viscous 
matrices constructed in such a way, the a and f n i vectors in 2.18 have to be ordered 
with the + and — modes alternating, i.e. a = {a*, <*7, a t > • • •}• 
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4S+1 
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a - 3- 3 
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a -30- 
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a ~32 
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a >4-l 
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a +41 

*♦4 3 

a +42 

a +42 *' 
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a ”4 1 

a -42 
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0 ' 
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a +5-2 

a + 5-2 

a : 5 -3 

a + 5-l' 

••.-iso- 

***.a+ 50 
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a +51 

a : 52 

a +52*'' 
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+ 

a -5-l 
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a“ 5 i 

a -52 

a ~52 

0 
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Figure C.l: Mass matrix as defined in 2.18. It is assumed that a second degree spline 
has been used, but the dimensions are shown for a spline of degree S. 


C.l. 2 Viscous Matrices 


In a similar way to the mass matrix, the viscous matrix is assembled from the follow- 
ing: 

•k z ^ 0 , kg 7 ^ 0 : 


K 

b: 



k 2 ((kg 2 - l)m 3 + m 4 ) + (kg 2 - l) 2 mi + (2kg 2 -f l)m 2 + mi 0 + 


(1 - kg 2 ) bt, - (2 + & 0 2 )bt 2 


— -=r- j2/c 2 4 m 3 ( — 3 + 4 kg + 2kg 2 — 4 kg^ + & 0 4 )kV 7 -(- 

3 k z 2 ((kg — l) 2 rri] 4- m 2 ) + (3 — 4 kg + 2kg 2 )tn% 4- mg + 

k z 2 ((l — Arg)bt] — 2bt 2 ) + (— 2 + 3 kg — &£ 3 )bt 3 + 

(kg - l)(bt4 - (kg - l)bt 5 )| 
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B+ = - 77 - \k z 2 ((k e + l)m 3 + m 5 ) + (k e - l) 2 ((ke + l)mi + mu) + 

Re 1 

(k$ + 2)m2 + mi3 — (2 + k$)ht2 
1 2 

B+ = — k z 2 ((ke + l)m 3 + me) + (ke — 1 ) 2 ((k$ + l)nii + 11112 ) + 

(ke + 2 )m 2 + m 14 - (1 + fc*)bt 2 

• k z = 0 , k Q ^ 0 : 

B+ = — ~~ ^ 2mi + l ) m 2 + m io + (1 "■ A:^ 2 )bti — (2 + k$ )bt 2 

BI = — ^(—3 + 4 k$ + 2^ 2 — ike 3 + fc^ 4 )ni7 + (3 — 4 ke + 2^ 2 )ni8 + + 

(—2 + 3 ke — &0 3 )bt 3 4- (k$ — l)(bt 4 (k$ 1 )bts) J 

Bl = B~ = 0 

%k z ^ 0, ke = 0: 

k 2 r 1 

B+ = Bl = B; = --^[^ 2 m 3 + mi +m 2 -bt 2 j 

BI = -4-[2^ 4 m 3 - 3(m 7 - m 8 ) + 3^ 2 2 (mi + m 2 ) + m 9 + 

Re L 

^ 2 2 (bt! — 2 bt 2 ) — 2 bt 3 — bt 4 — bt 5 ] 

•k z — 0 , kg = 0 : 

B+ = + m 2 - bt 2 

BI = — 2- fm 8 — m 7 — bt 5 — bt 3 

Re L 

Bl = b; = 0 
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C.2 Nonlinear Term 


Computing the nonlinear term is by far the most computationally intensive operation 
of the code; typically it requires over 85% of the computing time. This is due mostly 
to the radial direction which is a “slow” direction (requiring N? operations), whereas 
the azimuthal and streamwise directions are “fast”. 

The first step is to compute the vorticity and velocity components. Using the 


vectors in table B.l, 

(r,6,z,t) = Ai g,(r) + B t rg/(r) + C t (C.l) 

r L r 

= E l g l {r) + F l + G ig ,''(r) (C.2) 

u,(r,0,z,t) = H l ^ + (I l + 2J l )g, , (r) + J l rg l "(r) (C.3) 

v r (r,6,z,t) = Ktg,(r) (C.4) 

v e (r,8,z,t) = Jirg/(r) + Iig t (r) (C.5) 

v z (r,6,z,t) = -Gig t '(r) - Fi^^- (C.6) 

r 

where, the summation convention was used for repeated indices (a/6/ = a/6/), and 

the following were defined 


Ai(6,z,t) = EE-*V(a, + +ar)c*‘ (M+Jt '* ) 

kg k z 

dropping the summation signs and the exponential 


B t (9, z, t) = 

( —ik z 2 af 


\ 0 ,k 6 

Ci(0,z,t) = 

-ik e (k s - l)af 

D[(6,z,t) = 

ikgaf 

II 

uT 

5 

k z 2 (k e af + aj ) 
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Fi{0,z,t) 


Gi(6, z,t) 




h(0,z,t) 


(kg - 1 )a t 

| -«r 

[ 0 , kg = k z = 0 

f & 2 (1 - fc*)((l + + af) 

| (1 - kg 2 )o/ , k z = 0 

| k z (af + af) 

I a, + , k z = 0 


Jl(0,z,t) 


k z o/ , kz / 0 and fcj ^ 0 
a/ - , & 2 = 0 and ^ 0 

0 ,£:# = () 




—ik/kgaf + a, - ) 

—ikgaj , fc 2 = 0 and fcg ^ 0 

0 , fc 2 = kg = 0 


Once all those terms have been computed, they are transformed to physical space 
where the second step consists in computing the product W//* • (v x w). Defining the 
following three-dimensional (nonlinear) matrices 


fl = 

rR2 

1 9,,9,9 m dr 

& 

II 

* 

to 

25 

<cT 

3 

9- 

f3 = 

rR2 

j 9,/9,9 m dr 

II 

rR 2 

j g/g/dmdr 

tR2 

fs = j o 9,, 9,9m' dr 

<5T» 

II 

rR2 

1 9, , 9 ! 9 m' dr 

f 7 = 

rR2 

/ 9 / 9 , 9 / dr 

Jo 

rR2 

fs = / 9/ 9, 9 m ' dr 

J 0 

f 9 = 

rR2 

/ 9, ,9, 9 m" dr 

Jo 

II 

o 

cr 

rR2 

1 9,i Si' 9 m " dr 

rR 2 

fn =J o 9 /9,9 m" dr 




then, we also define: 


7,t = (E^Fl + Hmh) fl + [EraGl + H m Ji) rf 2 + 

( Im + 2 Jm) (I i rf 5 + Ji r 2 f 6 ) + F m ^F[ ^ ^ + Gi ^f 6 - + 
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G m {F, f 9 + Gi rf 10 ) + J m (I t r% + Ji r 3 f 10 ) 

0v = {A m h — E m Ki) f, 4- (C m Ii + F m Ki) 4- A m J, rf 2 -f C m J\ — + 

r 

D m Ji f 6 - G m l<i f 9 + Bm(h rf 5 + J, r 2 f 6 ) + (D m I, - F m I<,) - 

r 

+ Gi — ^ + B m Gi r 2 fg + 
(B m F + ImKi + 2 J m Ki) rf 5 + £> m ^ + Gi f 6 j + r 2 f 9 

f ,, = (v4 m f) + H m Ki) rf 3 + A m Gir 2 f 4 + C m (f, ^ + G, f 4 j + F m G, r 3 f 8 + 

(B m Fi + I m Ki + 2J m /\/) r 2 f 7 + Z) m (F/ f 7 H- G/ rf§) -f J m Ki r 3 f u 

S v = (A m Ii~E m K l )r { 3 + (C m I l AF m K l )- + C m JiU + A m Jir 2 UA 

V 

{D m 1 1 — F m Ki)i 7 + B m {I i r 2 t 7 + J[ r 3 f 8 ) + D m Jt rf 8 — G m Ki rfn 

where all the products f, r n are done under the integral sign, and the m and / indices 
have been summed; for example 

Nr Nr r R 2 

9„9,V rUr (C.7) 

m l Jo 

This means there are in fact twenty four different integrals, which are also pre- 
computed at the beginning of the code (some of them can be simply obtained by 
symmetry). Once 7 ; t, , Tj7, T;» and 8 /» are computed, they are transformed to 

Fourier space ( 7,t = 7 ( t ( kg , k z ), . . .) where the final answer is obtained: 

fnU' = 7T~r [ f f W;/ + • (v x Uj)e~ l ^ ke ‘ 8+kz ' z) rdrdOdz 

ZtcL 2 Jo Jo Jo 

= I + (ik el t (C.8) 

fnU 1 = F~7~ f f f w i'~ ' ( v X w)e~ l(ke ' 9 + k ‘' z) r dr d.9 dz 

ZttL z Jo Jo Jo 


Tp = (AmF + H m Ki ) fi + A m G,r f 2 + C m [F 4 
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= I-(i 7l t - T,7 ) + F- ft - G~ Sr 


- 0 and kg ^ 0 
= 0 

, kg = k z = 0 

C.3 FFT’s and Modal Reduction 


where 


I + = 


kz 

1 , k z — 0 


J + = 


I = k z , F = kg - l , 


kz 

1 , 
0 , 

G" = 


{: 


(C.9) 


In section 2.3.4, the concept of modal reduction was introduced to alleviate the con- 
centration of modes near the origin. Implementing this procedure in the context of 
the computation of the nonlinear term requires special attention. Consider figure C.2, 
which consists of an equispaced grid with 18 b-splines. 



Figure C.2: Modal reduction and Fourier transform size; the ordinate represents the 
number of Fourier modes. 


Azimuthal modes which fall below the staircase shaped dotted line represent the 
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only non-zero modes. Computing the product C.7 requires all the terms to be trans- 
formed on a physical grid with an equal number of points. It is clear that transforming 
all 18 modes on a grid containing N 2 points, would result in a large waste of com- 
putations, since roughly half the values to be transformed are zero. To reduce the 
length of the FFT’s, radial positions are defined where the lengths of the transforms 
are cutoff: for example, in figure C.2 the radial position equivalent to spline 11 is the 
cutoff value. Any mode lower than 11 would be computed on a grid of Ah, and from 
11 onward, on a grid of N 2 . 

Even though this approach is simple in principle, special care must be taken when 
approaching the overlapping region near the cutoff position. As seen in section A all 
integrals of the form shown in C.7 have support on 25+ 1 intervals, which means that 
computing the nonlinear term for /' = 8 requires information from 5 < V < 11. This 
can be a problem since modes 5-10 lie on a different grid than 11. The trick is then 
to FFT modes 5-10 to Fourier space with Ni modes and back to physical space again 
but on a grid of size N 2 . Only when this is done can modes 8 and up be computed. 
For V < 7, all the computations are performed on a grid of N\ . 

In practice, the code allows for up to six different cutoff levels; more than six 
shows no significant gain due to the extra work required by the remapping near the 
overlapping regions. 


C.4 Time Advance 

When the time advance method given in 2.30-2.32 is applied to 2.18, it becomes 


Ai«' = Bio:,, + A<7if n i(a n ) (C.10) 

A 2 a" = B 2 a' + At [7 2 f n i(a') + Cifni (<*„)] (C.ll) 

A 3 a n+1 = B 3 « ,/ + Af[7 3 f nl K) + C2f„i(« / )] (C.12) 

where A, and B, are respectively the effective mass and viscous matrices defined as 
A,- = A — /?,A<B and B; = A + o;AfB (C.13) 
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To minimize storage, the above method is implemented in the following way: defin- 
ing « 0 id and a act to be two storage variables the following algorithm is implemented 

• From the previous time step, set a a ct = c*old = <*n- 

• First RK sub-step: 

1. Compute the new At 

2. Q! ac t = fnl(^act) 

3- QJold = B]QJ 0 ]d 

4. a ac t = a old + 7 i A<a act 

5. <*old = a V~^T A ' i (= fnl(«n)) 

6. Enforce the boundary conditions (BC) on « ac t 
I- cx ac t = Aj a a ct 

8. «old = B2« a ct + Cl AtO: 0 ld 

9. Enforce BC on ct 0 id 

• Second RK sub-step: 

1 ■ C^act r fnl(^act) 

2. a a ct = a 0 id + 7 2 Ata ac t 
3-«oid = ^^ (=U<*')) 

4. Enforce BC on a ac t 
5- C*act = A 2 <*act 

6. <* 0 ld = B3Q!act + C2A<Q:old 

7. Enforce BC on o: 0 id 

• Third RK sub-step: 

1- <^act = fnl(^act) 

2. Ctact - <*old + 73 Atttact 
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3. Enforce BC on at acl 

4- <^act — -A. 3 C*act ( = ®n+l) 

5- «old = «act (= Ct n+ 1 ) 

• Proceed to next time step 


Constant mass flow 

In order for the mass flow to be constant, a pressure gradient must be applied in 
th'e streamwise direction on the mean flow, i.e the (0,0) mode. Let us symbolically 
rewrite any of the steps in C.10- C.12 as 

£o*„ +l = RHS + (wj' - (r;0,0),/e z ) (C.14) 

where RHS refers to any right hand side, fe x is the streamwise pressure gradient 
and v z = v 2 (r, f;0,0). Since this pressure gradient is not explicitely known, Moser 
(Private communication) suggests the following: split the solution step such that 


= RHS (C.15) 

Cv Zn+1 = (w,,"(r;0,0),e z ) (C.16) 


where the e z on the right hand side of C.16 refers to a unit pressure gradient. Once 
Vzn+i an< ^ Vz n + 1 bave been computed, the actual v ZnJtl can be obtained by 


U *n+1 = ^n+1 + ^ 


z n + 1 


(C.17) 


All that remains to be computed is S ; defining 

rR 2 




/*« 2 

= / v z 
J 0 


r dr 


(C.18) 


which is no more than the integral in 2.3 that has to be made constant Vf, and 

r/? 2 


m 


[K 2 

"Jo 


~ rR 2 


r dr 


r dr and 


(C.19) 
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then, 

6 = C * ^- m (C.20) 

rh 

Finally, note that splitting the solution in two steps does not require more memory: 
since the imaginary part of the (0,0) mode is zero, we can use that memory location 
to solve for C.16. 




Appendix D 

Benchmark Solutions 


In this appendix we present benchmark problems which were used to validate the 
computational procedure: Stokes problem and linear stability theory. 

D.l Solution to Stokes Problem 

Stokes flow represents the low Reynolds number limit of the Navier-Stokes equations. 
Since diffusion dominates convection, it is possible to neglect the nonlinear convective 
term from the equations rendering them linear. Defining A* to be the eigenvalues 
(assuming the solution obeys u ~ e* J< u s such that d/dt = X s ) and u s and p s the 
eigenfunctions of the Stokes equations, 

V • u s = 0 (D.l) 

A s u s + S7 Ps = V • Vu s ‘ (D.2) 

net, 

with u s = 0onr = i ?2 = l* 


D.1.1 Exact Solution 

This problem was solved by Salwen and Grosch [40] by means of vector and scalar 
potentials. Their solution was restricted for k z ^ 0 but was extended to include all 
wave numbers. 
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• |k z | > 0 

The eigenvalues are 

where the are the roots of the following transcendental equation 

J J LW\ ( JLW , W \ kAh 1 + P 2 ) 
\JM) \BJtM k,h,(K)) *» 2 /? 2 ' 


(D.3) 


(D.4) 


and Jk„ are Bessel functions and I k$ are modified Bessel functions both of the first 
kind (see Abramowitz and Stegun [1] chapter 9). Also, it is implied in D.4 and what 
follows, that k z = | k z | since the equations do not depend on the sign of k z . Since the 
eigenvalues X s represent the decay rate of the eigenfunctions, they are all real and 
negative. The eigenfunctions were obtained using Mathematica and are given by 


kg , k z ) kzl kg {k z r) -(- Q,Jk e +i (/? s r) bJ kg —]({3 s r) 
h e (k z r) 


and 


Ug s (r-,k e ,k 2 ) = k B 


j kg , k z ) — k z I kg (k z r) 4“ 

k* 


~ aJ ke+1 (P s r ) - 6«7jt e _i(/3 i r) 


(fa + l) Jfa+l(Ar) +AJl. +1 (ft r ) 


(D.5) 

(D.6) 


+ 




where 


a — 


_ kjhJM- kJlJk,) 


W*,+i(A) 


, wuy+wyy 

2 Jt,-i(A) 


(0-7) 


( 0 . 8 ) 


(D.9) 


When kg = 0, D.6 yields a trivial solution (ug 3 = 0); for this case the following 
form is used 

Ug,(r]0,k z ) = Ji(0 t r) (D.10) 


where /? s = j 1)S and j n<s is the s th zero of the Bessel function of the first kind of degree 
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Tl (l.C. ) fb 1,2,...). 

•k z = 0 

The eigenvalues and eigenfunctions are 


\ Jjc g+l^s 

s ~ Re b 




(Dll) 

(D.12) 
(D. 13) 


When k z = 0 the streamwise component decouples from the other two components, 


u*,(r; 4,0) = Jk e (jk a ,.r) 


(D.14) 


with the eigenvalue A 5 = 


D.l. 2 Numerical Solution 


The numerical approximation to this problem is identical to what was developed in 
2.18-2.20. With the mass and viscous matrices computed and assembled as shown in 
appendix C, a generalized eigenvalue problem is obtained by replacing a by Aa; 2.18 
becomes 

AAar = Ba (D.15) 

where a is the eigenvector and A the eigenvalue. The eigenvectors are in fact b- 
splines coefficients that can be summed following 2.16 to yield the corresponding 
eigenfunction. D.15 was solved using MATLAB which implements the EISPACK 
libraries. 
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D.2 Linear Stability 


Linear stability theory concerns itself with equilibrium solutions to the Navier-Stokes 
equations and examines the effect of growth or decay of perturbations (see Canuto 
[7]). Let u and p be solutions to the Navier-Stokes equations; define U and P as 
the equilibrium solutions and u and p as the perturbations such that the velocity 
and pressure are given by U 4- u' and P + p' ■ By substituting these definitions in 
the Navier-Stokes equations (2. 1-2.2) and neglecting quadratic terms, the linearized 
Navier-Stokes equations are obtained 


V u = 0 (D.16) 

^ + U- Vu + u- VU + Vp = -2—V • Vu (D.17) 

ot R^b 

In studying stability of pipe flow, it is customary (see for example Salwen et al. 
[40, 39]) to adopt 

U = (1 — r 2 )e z (D.18) 

with a forcing term as defined in 2.2 given by / = 4 / Re^, i.e. this value of f ensures 
that the mean flow D.18 is preserved. By assuming that perturbations have a solution 
of the form 

u = u(r)e‘ (M+ *‘* ) - iu " (D.19) 

D.17 can be written in a form similar to Stokes’ problem 

Au + U • Vu + u • VU + Vp = — — V • Vu (D.20) 

Rcb 

where A = —iuj. Perturbations are said to be stable if Im(a;) < 0. Equations D.20 is 
now a generalized eigenvalue problem that is solved similarly to D.15 


AA« = (B — 0)c* 


( 0 . 21 ) 
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with the added matrix O defined as 

O = / w^-fU-Vu^ + u^-VUjrdr (D.22) 

Jo 

As for Stokes flow, the eigenvectors a are b-splines coefficients that can be summed 
following C.4-C.6 to yield the so-called Orr-Sommerfeld waves. 

Because D.22 represent a term which must be added to the computer code in order 
to solve D.21, it must also be validated. Leonard and Wray [20] computed A with 
their fully spectral code for Re b = 9600, k e = k z = 1; by ordering the eigenvalues such 
that Re(A!) > Re(A 2 ) > . . ., they obtained A! - -0.023170795764 -i 0.950481396668 
with N r = 37 Jacobi polynomials. This value compares well with Salwen et al. [39] 
who obtained Aj = —0.02317 — i 0.95048 by solving D.16 and D.20 using expansions 
in terms of Stokes eigenfunctions. 

Table D.l gives the difference between the b-splines and the Jacobi polynomials, 
using the same number of functions (albeit on an equispaced grid for the splines). It 
is interesting to note that for 5 = 7, the fully spectral Jacobi polynomial code and 
the b-splines method give virtually the same answer. 

Table D.l: Relative difference between b-splines of different degrees and Jacobi poly- 
nomials (N r = 37). 


5 

— Re(Ai) 

— Im(Ai) 

1 

3 

0.023169 

0.9504816 

1.51 x 10" 6 

4 

0.0231708 

0.95048142 

2.06 x 10" 8 

5 

0.0231707960 

0.9504813968 

3.35 x 10" 10 

6 

0.023170795752 

0.950481396672 

1.31 x 10- 11 

7 

0.023170795766 

0.950481396674 

6.75 x 10~ 12 
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